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Bemerkungen zur numerischen Integration 
gewohnlicher Differentialgleichungen n-ter Ordnung” 


Von 
HEINZ RUTISHAUSER 


Es ist bekannt, daB man eine Differentialgleichung n-ter Ordnung in ein 
System von ” Differentialgleichungen erster Ordnung verwandeln kann; davon 
macht man insbesondere fiir die numerische Integration gerne Gebrauch. Es 
besteht aber auch die Méglichkeit, die numerische Integration einer Differential- 
gleichung m-ter Ordnung ohne Verwandlung in ein System durchzufihren, nam- 
lich mit den Formeln von Kutta-Nystr6 [1], bzw. KUTTA-ZURMUHL [2] und 
ADAMS-FALKNER [3]. 

Es wird in diesem Bericht versucht, die beiden Méglichkeiten — direkte 
Behandlung einerseits und Verwandlung in ein System andererseits — in bezug 
auf ihre Genauigkeit zu vergleichen. Zwar wurde eine Untersuchung iiber den 
Fehler der Kutta-Zurmiihl-Formeln bereits von L. COLLATZ und R. ZURMUHL [45] 
angestellt, doch bezieht sich dieser nur auf einen einzelnen Integrationsschritt 
und nicht auf den Fehler im Gesamtverlauf der numerischen Integration. 

Der Vergleich der verschiedenen Verfahren erfolgt an Hand von linearen 
Differentialgleichungen mit konstanten Koeffizienten, fiir die man die durch die 
numerische Methode induzierte Differenzengleichung exakt lésen kann. 

Die durch Rundungsfehler verursachten St6rungen sind hier nicht beriick- 
sichtigt ; vielmehr wird angenommen, daB alle bei der numerischen Integration 
auftretenden Rechenoperationen exakt ausgefiihrt werden. Dies ist in Anbetracht 
der verhaltnismaBig groBen Stellenzahl der meisten Rechenautomaten durchaus 
realistisch. 


§ 1. Grundlagen der Fehleruntersuchung 
Wenn wir mit einem Integrationsverfahren vom Typus RUNGE-KUuTTA oder 
mit den Formeln von KutTa-NystR6OM bzw. KuTTA-ZURMUHL eine lineare Diffe- 
rentialgleichung -ter Ordnung mit konstanten Koeffizienten: 


yi) + cy yD + coy) 4.0 +o, sy’ +e, y =0 (1) 


numerisch integrieren, so sind die Werte von y, y!,..., y"- am Ende eines 
Integrationsschrittes linear von den entsprechenden GréBen am Anfang desselben 
abhangig, d.h. wenn man die  GréBen y, y!,..., y"-) zu einem Vektor (x) 





* Eine Zusammenfassung (ohne Beweise) dieser Arbeit ist in ZAMP 6 (1955), 
p. 497—498 erschienen. 
Numer. Math. Bd. 2 20 
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zusammenfaBt und in iiblicher Weise /, als Abkiirzung fiir den Wert von /(*) 
an der Stelle x,=%9 +, setzt, so gilt: 


Dats =C- 5, (2) 


wobei die Matrix C noch vom numerischen Verfahren, von der Differentialglei- 
chung und von der Schrittweite  abhangt. 


Es ist bekannt, daB die exakte Lésung der Gleichung (1) bei gegebenem 
Anfangsvektor ) (%)) =) durch 


(x) = e*—*4 pg 


gegeben ist, wobei 
0 1 0 0 
0 0 1 0 ‘ (3) 
A= .. 
0 0 0 4 
—Ch —~—Cy-y —Cy—g ae 





Insbesondere ist also y, =v (x,) =e*"4 9, so daB der Fehler nach k Integrations- 
schritten offenbar durch 


D, — 9 (%,) = {C* — (€*4)*} 9 (4) 


gegeben ist. Als MaB fiir den Fehler wahlen wir freilich nicht einfach diese 
Differenz, sondern eine GréBe, die im wesentlichen den relativen Fehler in y 
wiedergibt. Dabei wird zweckmaBig noch durch die Integrationsstrecke dividiert ; 
wir erhalten dann als Fehlermaf F: 


In y(¥) —In 9 (4) 
F= —z (5)? 





(relativer Fehler pro Langeneinheit der Integrationsstrecke). 


Das asymptotische Verhalten (fiir x->0o) der FehlergréBe F 14Bt sich in der 
Regel leicht berechnen: Gibt es einen Eigenwert A der Matrix C mit gréBtem 
Absolutbetrag und einen Eigenwert « von A mit gréBtem Realteil (oder eventuell 
je ein komplexes Paar AZ, «% mit diesen Eigenschaften), so gilt praktisch fiir 
jede Lésung y(x) von (1) und jede numerisch erhaltene Lésung y, bei groBem x 
bzw. groBem k: 


y(x)res wd. (6) 


Man kann beide GréBen unter Beachtung von x=kh vergleichen und findet 


Iny~rwax; Inywkinda= “Ind; 
somit (fiir groBes x). (7) 


ah—Ind 
Fw ” ar 





? Durch numerische Integration erhaltene Werte sind immer durch eine Tilde 
gekennzeichnet. 


® Bei oszillierenden Lésungen ist dieser Ansatz sinngema4B zu modifizieren. 
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§ 2. Die Fehlerordnung eines Verfahrens 


In allen praktisch interessanten Fallen ist es so, daB das (noch von x abhiangige) 
FehlermaB F bei festem x mit 4 gegen 0 geht, und zwar wie eine gewisse Potenz 
von h: 

Definition. Sei ® eine Klasse von Differentialgleichungen und V ein nume- 
risches Integrationsverfahren. Die gréBte Zahl p mit der Eigenschaft, daB fiir 
alle Differentialgleichungen aus & bei numerischer Integration mit V der Fehler 
y — y baw. (| — | bei Systemen) bei festem x und fiir 4-0 stets von der GréBen- 
ordnung O(h?) ist, nennen wir die Fehlerordnung des Verfahrens V. 

Wir wollen hier ausdriicklich festhalten, daB sich die Fehlerordnung in den 
nachstehenden Untersuchungen immer auf die Klasse der linearen’ Differential- 
gleichungen mit konstanten Koeffizienten bezieht. Damit sind die in dieser 
Arbeit abgeleiteten Fehlerordnungen obere Schranken fiir die Fehlerordnungen 
in bezug auf jede umfassendere Klasse von Differentialgleichungen. 


Hieriiber gilt : 
Satz 1. Hat die zum Verfahren V gehérige Matrix C fiir jede feste Diffe- 
rentialgleichung (1) die Eigenschaft, daB fiir ein gewisses p>0 


li C—ehA 
ee 


=K (8)! 


existiert, und ist fiir wenigstens eine Differentialgleichung K +0, dann hat das 
Verfahren V genau die Fehlerordnung . 


Beweis. Bezeichnen wir e’4 mit B, dann lauft die Behauptung offenbar 
darauf hinaus, daB C*— B* fir h->-0, hk=x=const die GréBenordnung O (h?) 
hat. Nun gilt allgemein; wenn man C — B=<¢ setzt: 

Ck— BE=C**24C*-*¢B+4C*-%.¢. Bt+.--4+CeB*+¢B, (9) 


Fir h—0, also = =k->oco gilt wegen (8) e—0 und sogar C/— B/-—0 solange 


7k. Damit wird die rechte Seite von (9) in erster Annaéherung zu 


x 
7 feo *€- elf —*0) A dé, 


%o 


und somit gilt nach (4) asymptotisch fiir kleines h: 
(x) — 9 (x) ~W( fel™P4 K et) 4 de) ny = WP 4 (2), (10) 


wobei der Vektor 3(x) die folgende Differentialgleichung erfiillt: 


; § (x) — A 4(x) = K(x) = v(x) 
mit (41) 
§(%) = 0. 


Man muB jetzt offenbar nur noch beweisen, daB die Lésung von (11) nicht 
fiir alle Anfangsvektoren bj, und fiir alle Differentialgleichungen (1) so beschaffen 





? Dabei ist A die in (3) definierte, zur Differentialgleichung gehérige Begleitmatrix. 
20* 
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sein kann, daB ihre erste Komponente z,(x) identisch verschwindet. In der Tat 
ist dies bereits fiir die Differentialgleichung y”) — y=0 unméglich, bei welcher 


04 0 
001 0 
A= at 
00 o41 
‘es 


In diesem Fall lautet namlich (11) in Komponenten geschrieben: 


, 

2 — 2 = 0, (x) 
’ 

Zg — 2 = V2 (x) 


2, — % =,,(%), 


wobei die v;(x) die Komponenten des Vektors »(x)=Kvy(x) sind. 


Unter den Zeilen der Matrix K sei nun die j-te die erste von Null verschiedene. 
Dann ist v,=v,= --- =v;_,;=0, aber bei geeigneter Wahl von jj ist v;(x)=0 
und insbesondere v,(%9)=+-0. Damit fiihrt die Annahme z,=0 auf Grund der 
Gleichungen (12) sofort zu einem Widerspruch, indem dann namlich wegen (12) 
2, =%, = --- =z; =0 sein miiBte, wahrenddem aus der j-ten der Gleichungen (12): 
2; —2;4,=0;(x) wegen 2;,1(%») =0; v,;(%») +0 sofort z;(x%») + 0 folgt, wzbw. (Fiir 
j=n ist z;,, durch z, zu ersetzen.) 


§ 3. Numerische Integration mit der Eulerschen Formel 


Fiir ein System von Differentialgleichungen ’=/(x,) lautet das Eulersche 
Verfahren wie folgt: 


Beginnend mit jg) (%,) berechne man fiir k=0, 4, 2, ... 
Xp, 0,) > Dp 

ap ee (13) 
De + hoy > 0p41- 


Dank seiner Einfachheit kann das Eulersche Verfahren, das sonst kaum ver- 
wendet wird, als Modell fiir die weiteren Untersuchungen dienen: 


Aus (13) entnehmen wir sofort, daB die in (2), §4 definierte Matrix C fiir 
ein lineares System )’= At mit konstanten Koeffizienten beim Eulerschen Ver- 
fahren wie folgt lautet: 


C=1+hA. (14) 


Da hier lim ee = — £ * offenbar existiert, hat das Eulersche Verfahren nach 
Satz 1 die Fehlerordnung 1. 

Will man eine Differentialgleichung m-ter Ordnung integrieren, so ist es nahe- 
liegend, das Eulersche Verfahren wie folgt zu modifizieren (in Analogie zu den 


Formeln von KuTTa-ZuRMUHL) : 
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Ausgehend von yp», yo, .-., y$"~» berechne man fiir k=0, 4, ... 


one Te Hk te => yf 


~tit 


Wthnt+ a ~ Fi iad = a > Vai 


~w 


t+ hyp eid i = ae => Vaya (45) 


sip + hoy => yf". 
‘ n—v h* ‘ es 
Allgemein: >’ _ yet = yf), (v=0,4,...,%—4). 


0 





Wendet man das modifizierte Verfahren auf die Gleichung (1) an, so hat die 
Matrix A die Form (3) und es wird 


h? he-} o 
‘> wem\ far 
1 a o 
C= 1 Ln? . ” 2 pay (Cys Cy—1, +++» Cg, C4). (16) 
Le cn y h2 
0 4h 2 
4 h 


Eine leichte Rechnung ergibt hier 


0 0 0 
— 0 0 aia 
e=C — eh => : : . + Gliederhéherer (17) 
0 0 0 Ordnung inh, 
— OC, Cy —%Cn-1> Co— cf 


so daB nach Satz1 auch das modifizierte Verfahren (15) die Fehlerordnung 4 
besitzt. Freilich ist hier v,=v,=---=v, -,=0, so daB man aus der ,,Fehler- 
gleichung*‘ (12) z,(%), 23(x), ..., z,(%) eliminieren kann; es resultiert dann folgende 
Differentialgleichung fiir die erste Komponente von z (x): 


2") +c, 2l™—D 4 cy 2-9) +... +, 4 2a + Cy 2 = Vy (%) 
mit (18) 
2 (%) = 2 (%) = +++ = 24") xy = 0. 


Gema&B der Bedeutung von z(x) folgt aus (18) sofort, naB fiir kleines x — x 
F(x) — y(x) whe, (x) wh EA). (oy — a) (19) 


Der Fehler bleibt also beim Integrationsbeginn zunichst noch sehr klein, so daB 
eine héhere Fehlerordnung vorgetauscht wird, wenn man nur wenige Integrations- 
schritte ausfiihrt. Wir wollen diese Erscheinung an Hand der Differentialglei- 
chung y'Y— y=0 mit y)= yo = ¥o = Yo. = 1 experimentell untersuchen. 
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Tabelle 1 * 
h=0,01 h=0,1 
x Exakte Lisung mit Formel (15) mit Formel (15) 
_—— EuLer (13) —-- EvLer (13) 
y=ex y | 10°F 10° F 7 10° F 10° F 
0 1 1 — 1 
0.2 1.221 403 1.221402 1 1.221400 9 
0.4 1.491 825 1.491 819 9 | 1.491777 80 
0.6 1.822119 1.822089 28 | 1.821851 244 
0.8} 2.225541 2.225441 56 2.224631 511 
1 2.718282 2.718026 94 4967 2.715916 871 46898 
1.5] 4.481689 4.480224 | 218 | (konstant) 4.467878 | 2058 | (konstant) 
2 7.389056 7.383 767 358 7.338740 3416 
3 20,085 54 20,04909 | 605 19,73654 | 5843 
4 54,59815 54,43012 771 52,989 62 7476 
5 148,41316 | 147,768 i4 871 | 142,25786 8472 
10 22026,47 21 794,34 | 1059 | 19865,18 | 10328 
15 3269 017 3214468 | 1122 | 2774178 | 10942 
| 1 
fo) 1247 | 4967 12172 46898 





* Samtliche Werte wurden mit 11 wesentlichen Stellen mit der ERMETH berechnet 
und nachher gerundet. 


Die Tabelle zeigt folgende Eigenschaften des modifizierten Eulerschen Ver- 
fahrens (15): 

(a Die Uberlegenheit gegeniiber den urspriinglichen Eulerschen Verfahren (13) 
fiir kleines x. 


b) Das starke Anwachsen der FehlergréBe F. 
c) DaB F ungefahr zu h proportional ist, d.h. daB die Fehlerordnung 1 ist. 


§ 4. Numerische Integration mit den Formeln von RUNGE-KUTTA 


Fiir ein lineares System y’= Av mit konstanter Koeffizientenmatrix A ergibt 
die Integrationsvorschrift von RUNGE-KUTTA!: 


i), = Ausgangswert fiir den Schritt von x, nach x,4,; =A, 
~ ~ h ~ h - Hod | ro 
by =e t+ 5 hk =(E +5 A) by = 4 ds 
a a. h ee - . 
Dax =n +> Oe =(E+5 A+ 4 A?)i, Yee =A Des 
~ ~ ~9 h? 2 h 3\¢ ~? ~ 
Deve =O, thigy =(E +A + > 4 i 4 A?) i, Dees = 4 Dawe 


~ ~ h  ~ ~? ~? ~? 

Det si De a 5 (0, + 204+ 20eet Sees) 
alE+h4+ 20+ 24 = As) § 
- a". 2 








1 Vgl. etwa [4], S. 33 oder [7], S. 344—347. 
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Somit ist fiir dieses Verfahren 





m me gay yo he ga 
C=E+hA+ = A+ 7 AS4+ 27 A (24) 
und damit 
- C—ehA As 
i (22) 


so daB das Verfahren nach Satz 1 die Fehlerordnung 4 hat. Man kann ferner den 
asymptotischen Wert der FehlergréBe F fiir x—> co bestimmen: 


Hat die Matrix A einen Eigenwert « mit gréBtem Realteil, so ist 
h? hs h* 
= ~_ es a oe ol 
A=1+ha+ 7 a ye 


ein Eigenwert von C, der bei hinreichend kleinem h iiber alle andern Eigenwerte 
dominiert, so daB sich dann fiir groBes x jede Lésung ) von y’=Av und das 
zugehorige Aquivalent f wie folgt verhalten (y steht fiir eine beliebige Kom- 
ponente von bh): 


y (x) ~e%* 
. h? h3 ht. \k (23) 
jp~(tthat s+ Fat+ at] , 


somit fiir groBes x und kleines h: 


2 3 4 


h oar 





Fw (24) 


Es sei aber ausdriicklich festgehalten, daB fiir zu groBes h méglicherweise 
der Fall eintreten kann, daB fiir einen zweiten Eigenwert B von A 


zwar Re (B) < Re(a), 
| (25) 


aber [1+hB +4 p*-+~ p+ ps >that Sort Sass H ac! 








In diesem Fall wird sich die numerische Lésung fiir groBes x wie 
h2 h® h4 k 
(1448+ 5 +S B+, BF) 


verhalten und daher stark von e** abweichen (Instabilitat). 


§ 5. Numerische Integration mit den Formeln von ZURMUHL 


R. ZURMUHL hat in [2] allgemeine Formeln zur numerischen Integration einer 
Differentialgleichung m-ter Ordnung (ohne Verwandlung in ein System von 
Differentialgleichungen erster Ordnung) aufgestellt. Diese Formeln enthalten 
als Spezialfalle fiir »=1 die Formeln von RuncE-Kutta (s. § 4) und fiir »=2 
die Formeln von Nystr6m [1] (s. auch § 6). 
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Wir wollen hier zeigen, daB auch das Verfahren von ZURMUHL fiir alle n 
die Fehlerordnung 4 besitzt, entgegen anderen Behauptungen, die den Zurmiihl- 
schen Formeln eine héhere Genauigkeitsordnung zuschreiben. 

Zunachst beschreiben wir die Vorschrift von R. ZURMUHL zur Integration 
der Differentialgleichung y” =/(x, y, y’,..., y"—"); der Schritt von x, nach %,4, 
lautet dabei folgendermaBen: 


Ausgehend von 9, 9,9; --. Vf"—") berechne man: 








a) F(X Ins Deo Th ver Ie) => yf" 
b) fiir »=0,1,2,...,.%—1: oa" =» 5) 
h ~ n— ay 
c) txt 51 Fe Fer TK) oH 
d) fiir y=0,1,...,%—2: ye) => yl), 
y=n—1: ye d+ — 2 5 »yoo" 
h ~ n— n 
°) 1+ Sam S85) 958 
f) fir y=0,41,....%—1: (26) 
> h* yet) 4 yn =» 
=, rm on HI Oe ae 
g) t (ptr Vawe . an »* re yee) > yn) 
h) fir v=0,1,...,m—1: 


(n—v)?® 54" +-2(0—2) P+ ILI +(2+9—M) Hoe _, po) 
(n—v+1) (n—v+ 2) 
n—v—1 a . 
2. = et +; = Rt) =>), 


x=0 











Es ist klar, daB die Fehleruntersuchung fiir diese Formeln wesentlich kom- 
plizierter ist als bei den Eulerschen Formeln und daher nicht mehr in der gleichen 
Weise wie in §3, d.h. an einer beliebigen linearen Differentialgleichung u-ter 
Ordnung mit konstanten Koeffizienten durchgefiihrt werden kann. Wir wollen 
uns vielmehr auf die spezielle Differentialgleichung y") —y=0 beschranken, 
welche aber bereits zeigt, daB die Fehlerordnung des Verfahrens von Kutta- 
ZURMUHL nicht gr6Ber als 4 sein kann. 

Wenden wir die Formeln (26) auf die Differentialgleichung y™ — y=0 an, 


so erhalten wir: 
nn 


= Y M—T poe an nn 2 oye ts 
* =Vae =>; — 5 Jaws = Dh + nt 2s Beet Hh 





und damit: 


(n+ 1 —») (n+2—0) RO = (n— 9)" + 4(n —») 5p + 
0 


+Q+9— md TIP +5 nl pat 2 x! 0} 
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Ordnet man noch nach Potenzen von A und beniitzt y= y"**), so erhalt man 
fiir n=>3: 

RW) —— a;(%) h ,)(%+1) h® (+2) 

=I + pian e+ epi tae + 








4+ =e ~(n eo ° ‘ 
+ he 12(n “4 3 a ») yet) + hdhere Glieder in h 


somit 
a+3— 5, (4+0— n) 
=) 0 yore) 4 Waistcoat? y+) 4+ hdhere Glieder. (27) 


Die Differenz zwischen diesem Ausdruck und der, ,exakten Formel‘‘y{’) ,= => “ft ) 
betragt somit in erster Annaherung 


pes 4+v—n oe 1 n+3—» O(m+3) __ ___(n—») ( (s—t—9) n+3—v5(n+3)_ 
R= \42(n+2—9)! (n+3—»)! )h Wk 12(n+3—»)! | h yy 





Fiir y=0, 1, ..., # —3 ist dieses Restglied von der GréBenordnung </*, aber 


5 
fiir yn —2 wird R,_.=— a y't®), wahrenddem R,_, wegen dem Faktor 


n —1—v verschwindet. Fiir y= —1 tragt also erst das nachste Glied der Ent- 
wicklung zum Fehler bei und dieses hat ebenfalls die GréBenordnung hi’. 

Damit ist gezeigt, daB die Differenz C —e*4 fiir die Differentialgleichung 
y™ — y=0 genau die GréBenordnung h5 hat, so daB das Verfahren von KutTA- 
ZURMUHL (fiir »=3) nach Satz 1 héchstens die Fehlerordnung 4 haben kann. 
DaB diese wirklich 4 ist (auch fiir »=2) geht aus den Ausfiihrungen bei COLLATZ 
[4], Kap. I, § 33 hervor. 

Genauere Berechnungen zeigen, daB die in Satz 1 definierte Matrix K fiir die 
Differentialgleichung y” — y=0 (mit =4) die folgende Gestalt hat: 


Fiir n=1 (RUNGE-KuTTA) = (— 199). 
Fiir »=2 (KuTTA-NysTROM) K = x_( " = a , 
480 
0 
Fiir n=3 (KuTTA-ZURMUHL) K= “in . 
c: 0 
0 0 
a " te > ° 0 0 
Fiir »=4 (KUTTA-ZURMUHL) K= 0 0 all 
‘les0 = («OO 0 0 


Dieselben Nicht-Nullelemente — 4/99 und 4/sgg9 treten auch fiir alle hdheren m auf. 

Die Kleinheit der Elemente von K fiir => 4 im Vergleich zu »=1 und n»=2 
legt die Vermutung nahe, daB das Verfahren von KuTTA-ZURMUHL allgemein 
etwas genauere Resultate liefere als das urspriingliche Runge-Kutta-Verfahren. 
Numerische Versuche zeigen jedoch, daB diese Vermutung nicht allgemein richtig 
sein kann (s. § 6). 





272 HEINz RUTISHAUSER: 


§ 6. Numerische Experimente mit den Formeln von NYSTROM und ZURMUHL 


Sei y’’—(a+)y’+aBy eine lineare Differentialgleichung mit konstanten 
Koeffizienten. Im Falle «+f sind e** und e’* zwei unabhangige Lésungen. 
Intégrieren wir nun diese Differentialgleichung numerisch nach den Formeln 
von Nystr6M}!, so erhalten wir nach einiger Rechnung die Matrix C, die den 
Ubergang von §,=(,, 94) 2U Hyir=(Veir, Vers) Vermittelt, wie folgt: 


Ks Ses 
(en ci). 
wobei 
Cr =1— ap —™ (tp + opt) —% (eB + 02 Bt + a BY) +e (08 Bt + a2 


Cy =hte steeleesians 
+ Fy (0 + 08f + af? + BY) — Te (a + 2086" + a) 
Cxy= —hap—™ (02 B +0 8%) —* (a8 B+ 02 B+ 0.68) — 2 (a p-+08 f+ 02 B+ 0B) + 


+o ( (3.048? + 50383 + 3.0288) — Js Eisbaed 

Coz = 1+ h(a +B) + * (a? +ap +B) + Flaine 
4 © (at + ‘gs dite ith © (atB + 2086? + 2026? + aft) + 
+ ain (a4 B? +- 203 B% + ao? BS). 


Wir kénnen nun das asymptotische Verhalten der FehlergréBe F gemaB § 1, 
(7) aus den Eigenwerten der Matrix C bestimmen: Ist Re(«)>Re(), so wird 
sich die exakte Lésung fiir groBes x wie e** verhalten (von instabilen Lésungen 
abgesehen), so daB der absolut gréBere Eigenwert von C fiir kleines h in der 
Nahe von e** liegen muB. 

Wir kénnen ihn mit Hilfe des Newtonschen Verfahrens genauer berechnen; 
man erhalt offenbar 


D (eha) 
a D’ (eh) 


Cyi—A Cie 


ha 
Awe CG. Ca-i 


mit D(A) = 
Eine kurze Rechnung liefert 


D(e*) = a (4a? + 1108 -+ 106%) +47... D'(e) = h(a —B) 42%... 


Damit ergeben sich dann A und nach (7) der asymptotische Wert der Fehler- 
gréBe F (fiir groBes x und kleines h): 


h* at _ 408+ 110 B+ 106? 
480 a—B 





Fi & (28) 





1 Diese sind identisch mit dem Spezialfall »=2 der Zurmiihlschen Formeln (26); 
s. auch [7], S. 366 Formeln (4). 
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Dieses Ergebnis steht in Einklang mit der Behauptung, daB das Verfahren 
von Kutta-Nystr6m die Ordnung 4 hat. Fiir «—— erhalt man beispielsweise 
h* «5 
320° 
iiberlegen. Man beachte aber die gefahrliche Verschlechterung der Nystrémschen 
Formeln, wenn sich B und a nahern. Diese Verschlechterung wird durch numeri- 
sche Integration von y’’—2y’+y=0 mit y(0)=0; y’(0)=1 bestitigt (siehe 
Tabelle 2)?. 

Die fiir groBes x auffallende Verschlechterung der Nystrémschen Formeln 
veranlaBt uns, auch die allgemeineren Zurmiihlschen Formeln (26) in der gleichen 
Weise nachzupriifen. Zu diesem Zweck integrieren wir die Gleichung y!¥ — 4 y’” + 
6y"—4y’+y=0 mit y(0)=y’(0)=y"'(0)=0; y’”’(0) =6 mit den Zurmiihlschen 
Formeln mit »=4 (s. Tabelle 3?) 

Man beobachtet, daB F bereits bei x=10 etwa 100mal so groB ist wie bei 
den Formeln von RuNGE-KutTTA; die gefahrliche Situation bildet sich also bei 
den Zurmiihlschen Formeln noch wesentlich rascher aus als bei den Formeln 
von NystrOm. Gleichzeitig zeigt obige Tabelle aber, daB die Fehlerordnung des 
Verfahrens von KUTTA-ZURMUHL vermutlich auch in diesem kritischen Falle 
noch’4 sein diirfte?. 

Wir wollen indessen auch noch einige Falle betrachten, wo die Formeln von 
ZURMUHL wesentlich genauer sind, als die urspriinglichen Formeln von RUNGE- 
Kutta. Auf Grund der Ergebnisse von § 5 wird man vermuten, daB die Formeln 
von KutTTa-ZURMUHL tatsdchlich genauer sind (aber immer mit der Fehler- 
ordnung 4), wenn z.B. die Ableitungen y” und y’” nicht explizite in der Diffe- 
rentialgleichung auftreten. Das folgende Beispiel zeigt, wie man davon Gebrauch 
machen kénnte: 


| a) das Verfahren ist also in diesem Fall dem Runge-Kutta-Verfahren 


Es sei die Gleichung y”’ —1,6y’+0,6y=0 mit yy=yo=1 numerisch zu inte- 
grieren (exakte Lésung e*). Mit den Formeln von Nystr6m erhalt man bei 
Integration mit h=0,25:¥ (10) = 21 987,2 an Stelle des exakten Wertes 22026,47... 

Durch zweimalige Differentiation erhalt man aus der gegebenen Gleichung 
eine Gleichung 4. Ordnung, aus der man noch y” und y’” eliminieren kann. 
Man erhialt: 


mr 


y¥ —2176y' +1176y=0 mit yy=Y=%=Vo =1. 


Integriert man diese Differentialgleichung mit den Formeln von ZURMUHL eben- 
falls mit A=0,25, so resultiert y (10) = 22028,8. 


§ 7. Numerische Integration nach ADAMS-FALKNER 


Fiir jede hinreichend oft stetig differenzierbare Funktion y(x) gilt mit der 
iiblichen Schreibweise /, fiir /(%)»+ kh): 


N-1 
Ve=Ve-rth 2 Bx V" (v%) + Restglied, (29) 





1 Die Berechnungen wurden mit 23 wesentlichen Dualstellen auf der Z 4 durch- 
gefiihrt. 
2 Dies folgt tatsichlich aus den Ausfiihrungen bei CoLLatz [4], Kap. I, § 33. 
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wobei V” (y,) die v-te Riickwartsdifferenz der Funktion y’ an der Stelle x, bedeutet 
und das Restglied die GréBenordnung h" + hat. Die Koeffizienten B* sind durch 


die erzeugende Funktion 
co 


Le ** saa (30) 


0 


eindeutig bestimmt und bei CoLtatz [4], Kap. I, §15 angegeben. 

Auf der Formel (29), wobei N beliebig wahlbar ist, beruht ein Integrations- 
verfahren fiir Differentialgleichungen 1. Ordnung, welches die Fehlerordnung N 
besitzt (sog. Interpolationsverfahren von ADAMS-BASHFORT). Es kann sinngemaB 
auch auf Systeme von Differentialgleichungen und damit auch auf Differential- 
gleichungen m-ter Ordnung iibertragen werden. Jedoch muB man dabei fiir alle 
unbekannten Funktionen des Systems, bzw. fiir alle Ableitungen y’y’”’ ... y™ je 
ein Differenzenschema aufstellen. Demgegeniiber wurden von FALKNER [3] 
Formeln zur numerischen Integration einer Differentialgleichung »-ter Ordnung 
angegeben, bei denen man nur die Differenzen der héchsten Ableitung y” be- 
ndtigt: Die zugehérigen Interpolationsformeln lauten: 


n—o-1 N-1 
he f its 
y= >) = viet + ar-2 Dd Bs_ 00” (yf) +R, (Q=0,1,...m—1), (34) 
o=0 v=0 
wobei die Koeffizienten B%, bei CoLLatz [4] I, §16 definiert und angegeben 
sind. Die erzeugenden Funktionen sind hier 


~ 1 
p> Bu» *” = ese i P,( — In (1— x) (u =1,2,...,M), (30a) 


wobei das Polynom P,(x) vom Grade «—1 dadurch festgelegt ist, daB die nega- 
tiven Potenzen von x auf der rechten Seite von (30a) verschwinden miissen. 

Das Restglied R, der Formel (31) hat die GréBenordnung h"**~e; fiir e=0, 
also fiir die Berechnung des neuen y-Wertes, scheint der Fehler also von der 
GréBenordnung h"*+% zu sein. Wir wollen aber zeigen, daB die Fehlerordnung 
des durch (31) definierten Verfahrens dennoch nur N betragt, und zwar unab- 
hangig von n. 

Dabei ist natiirlich zu beachten, daB fiir den Integrationsbeginn erst einmal 
die Werte yy’... y"—- an den Stellen x %, ... Xy_» irgendwie beschafft werden 
miissen! und daB der ganze Integrationsverlauf und damit auch die Fehler- 
gréBe F davon abhangt, wie man dieses Anfangsstiick berechnet. Wir beschranken 
uns daher auf die Betrachtung des asymptotischen Verlaufs der FehlergréBe 
den man bei kleinem / nach den Methoden von §1 ermitteln kann. Sei 


L(y) = y + ay" +--+ a, 19’ +a, y =0 (32) 


die zu integrierende Differentialgleichung und es sei « diejenige Wurzel der 
Gleichung 
1 (ce) =a" + aya"? + ++ +a, 14 +a, =0 (33) 


1 Vgl. etwa CoLiatz [4], Kap. I, § 52, § 53. 
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mit dem gréBten Realteil und zwar sei dies eine einfache Wurzei. Dann verhalt 
sich fiir groBes x praktisch jede Lésung von (32) wie C -e** und fiir hinreichend 
kleines A verhalt sich dann auch jede mit Hilfe der Formeln (31) numerisch 
bestimmte Lésung y, fiir groBes k wie y A*, wobei A in der Nahe von e** liegt. 
(Das Verfahren ist namlich nach [6] fiir hinreichend kleines h stabil.) 


Mit dem Ansatz y,=y /*; y\?’=y, A* findet man durch Einsetzen in (31): 


n—o—1 


a N-1 v 
YM = ba a Yo+e r+ 2, Bro,» (1 die 7) Yn a 


o=0 


(34) 
(fiir 9 =0,1,...,%—1), 


denn die Riickwartsdifferenzenbildung entspricht bei einer Zahlfolge y A* einfach 
einer Multiplikation mit 1—>. Dividiert man (34) durch /*~?, so erhalt man 


die Gleichungen: 


n—o—l1 he N-1 . 1\" 
—Ant 2 at Vere + AW 5 BE-oe(t— 3) aie (35) 
(fiir 90 =0,1,...,%— 1). 
AuBerdem folgt aus (32) sofort 
Yn t+ UMYn—-1 t+ 42 Yuta ts + Anti V1 t+ an Vo = 9, (36) 


so da8 wir zusammen mit (35) +1 homogene Gleichungen mit »+ 1 Unbekann- 
ten Yo, 71, ---, Y, haben, die eine nichttriviale Lésung haben miissen. Dies fiihrt 
zur Determinantengleichung 


h2 pn-l AS 1\” 
yn-2 fics N-1 4\" 
—A+1 h Goa 47 Di Bt-a,(1— i) 
Dx(A)= —A+1 =0. (37) 
0 7 : N-1 ne 
—A+it Ah r,4——> 
2 Bi, ( i) 
a, ay-1 a,—2 a 1 








Da diese Determinante auf Grund der Definition der Bf, fiir A=e** exakt 
verschwinden mu8, wenn die Summen in der letzten Kolonne ins unendliche 
erstreckt werden (und dann noch konvergieren), kann man bei endlichem N 
und kleinem h die in der Nahe von e** liegende Wurzel A wie folgt angenahert 
bestimmen: 


Durch Entwicklung nach der letzten Kolonne bestimmt man erst Dy (e**); 
dann Dy (e**) ; alsdann ist angendhert (fiir kleines /): 


Pu). 


A ee et* — Di, (exh) (38) 
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Nun ist allgemein (falls die unendlichen Summen konvergieren): 


Dy (A) = Dy (A) + Dyan? h > pt (1 fide =) 


~ Dy ae Sp) 


Fa 


Dy. WS ft, (1 —3). 


(39) 





Dabei bedeutet D,, die durch Streichen der p-ten Zeile und q-ten Spalte aus 
Dy (A) erhaltene n-reihige Unterdeterminante. Es ist allerdings zu beachten, daB 
hier die Numerierung der Zeilen und Spalten mit 0 beginnt. 


Fir A=e*" wird D,,(A)=0 und 1— < hat selbst die GréBenordnung h. 
Da dann auBerdem sadmtliche Unterdeterminanten D, ,, von Dy (A) héchstens die 


GréBenordnung h”~1! haben, ist in (39) in erster Annaherung nur das Glied 


co Y D 
|. ae h > Bt, (4 — =) zu beriicksichtigen, sofern nicht jim — 1" —0. Nun 
N a od 


An-1— 
ist aber: 
h2 yn-1 
—A+1 h me FREER AEE (mat 
—A+1 h 

Dy-1,.0= on A f- 1 é ’ 

a+ h 

a, An-1 ae ay 








und dies wird fiir A=e*" in erster Annaherung, d.h. bis auf Glieder héherer 
Ordnung in h, zu (—1)"~1 h"-1 . a, a"~"=(— 1)" h"—! x", daa eine Wurzel von 
(33) ist. Damit wird nach (39) endlich: 
Dy (e*") ~ Bi (a h)**" (— 1)". (40) 
Nun folgt die Berechnung von —e. fiir A=": 


Wir erhalten durch Differentiation von (37) offenbar: 


dD 
a = — (Dog + Diy +++» + Dy-a,n—1) + 


+> (—° Doo $3 wad pe, (4 — i) |}- (41) 


Nun werden wir aber feststellen, daB alle Unterdeterminanten D,, wie h"~} 
klein werden, wahrenddem die D,_, , in der Summe rechts bereits abgesehen 
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von hk? schon den Faktor h"~1 enthalten, so daB diese Summe nicht mehr in Be- 
N-1 


tracht fallt. Wegen Ah’ B,,,(1—4;) ~2 wird allgemein 
0 











—ha h 
—ha 
0 A - 
om ha  ¥ 
Dg Aw | cnsonvssssessonnseesonnnessrsnssenssssaliestinannessectnsessesnaguscosensooeesnneesse 
—ha h 
—ha (42) 
0 0 Re 
: —hah 
By Anna Bn—0+1 | Op—o-1 a, 1 
is (— ha)? (— ert [s~e~4 + Kay, g-¢ + lle ao gprtns 
_— (— hy" [a"—* - oa + @,_-¢-1%"] 
und damit 


Doo t+ Dyr + +++ + Dy-1, n—1 = (— A)" [mat -? + (m — 1) aya? + oe + 
d 
+2ay_9-+ a1] = (— hy? <2 
Dabei ist y(«) das in (33) definierte Polynom. Somit nach (38), (40), (41): 
aN+a 


= — Big (43) 


BF. (ah)N+" (— 1)" 


ah 4 -* 
Am OO + (yet ye) 





und schlieBlich nach (7), §1 in erster Annaherung, d.h. fiir, kleines h: 


FE, = lim ~ BinaNt® pw 
um F(z) (a) h”. (44) 

Wir sehen damit: 

a) Die Fehlerordnung des Verfahrens (31) ist héchstens N. 

b) Fiir Differentialgleichungen (32), bei denen die Wurzel y(a«) mit gréBtem 
Realteil mehrfach ist, wird das Verfahren von ADAMS-FALKNER schlecht ab- 
schneiden (allerdings sind dann die obigen Fehlerbetrachtungen nicht mehr giiltig). 

c) Fiir eine Differentialgleichung wie y”’—y=0, bei der die Ableitungen 
y’ y”’,..., ¥"-» nicht auftreten, diirfte das Verfahren von ADAMS-FALKNER vor- 
teilhaft sein. Zum Beispiel ist fiir L(y)=y” — y offenbar x(a) =a" —1 und da- 
mit «=1 die dominante Wurzel. Ferner ist y'(a) =a"! und damit F,~ Fiz AN 
also gerade n-mal kleiner als beim gewohnlichen Verfahren von ADAMS. 


d) Die obigen Untersuchungen beziehen sich auf das Interpolationsverfahren 
von ADAMS-FALKNER, aber dieselben Erscheinungen zeigt auch das Extrapola- 
tionsverfahren. 
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Untersuchungen des Untergruppenverbandes 
endlicher Gruppen auf einer programmgesteuerten 
elektronischen Dualmaschine 


Von 


J. NEUBUSER 


1. Problemstellung 


Bei der Untersuchung gruppentheoretischer Probleme ist es oft niitzlich, die 
Fragestellung an Hand eines gréBeren Beispielmaterials zu testen. Die genauere 
Diskussion spezieller Gruppen bedingt jedoch, wenn es sich nicht um einen wohl- 
bekannten Typ, wie etwa abelsche Gruppen, handelt, einen unter Umstanden 
recht erheblichen Rechenaufwand. Es liegt daher nahe zu versuchen, elektronische 
Rechenmaschinen zur Untersuchung und Diskussion speziell gegebener Gruppen 
zu verwenden. 

In der Gruppentheorie liegen kaum systematische Methoden zur Struktur- 
untersuchung einer gegebenen Gruppe vor; haufiger benutzt werden nur einige 
Rechenverfahren wie etwa das der Restklassenabzahlung [1], S.12, fiir eine 
durch Erzeugende und definierende Relationen gegebene Gruppe, mit dem Ziel, 
eine Permutationsdarstellung fiir diese zu erhalten. Bei der Diskussion einer 
Gruppe macht man vielmehr gern Gebrauch von gruppentheoretischen Satzen, 
die gerade auf die spezielle Gruppe anwendbar sind, um langere Rechnungen 
méglichst zu vermeiden. Die Benutzung einer elektronischen Rechenmaschine 
wird dagegen erst dann sinnvoll, wenn man ein systematisches Verfahren fest- 
legt, welches auf beliebige Gruppen oder jedenfalls auf groBe Klassen von Gruppen 
anwendbar ist. In (durch die Speicherkapazitat der Maschine) begrenztem Um- 
fang kénnen dabei Anderungen der Rechenmethode auf Grund schon bekannter 
Teilergebnisse unter Benutzung gruppentheoretischer Satze durch Verzweigungen 
des Programms beriicksichtigt werden. 

Ein wichtiges Hilfsmittel, die Struktur einer endiichen Gruppe zu beschreiben, 
ist der Verband ihrer Untergruppen. Dieser kann durch die Eintragung der 
Klassen konjugierter Untergruppen und der Indizes der Untergruppen (voll- 
standiger Situationsplan) erganzt werden. Es ist zwar bekannt [3], daB es nicht- 
isomorphe endliche Gruppen & und § gibt, deren Untergruppenverbande V(@) 
und V() durch einen Verbandsisomorphismus aufeinander abgebildet werden 
kénnen, der Klassen konjugierter Untergruppen und Indizes erhalt; doch sind 
diese Beispiele fiir kleinere Ordnungen recht selten. Die Bestimmung des voll- 
standigen Situationsplanes ist um so mehr von Interesse, als man aus diesem 
auch wichtige charakteristische Untergruppen wie Kommutatorgruppe, Zentrum 
und Frattini-Gruppe ablesen kann. 
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Im folgenden soll ein System von Programmen beschrieben werden, das aus 
gegebenen Erzeugenden einer endlichen Gruppe nicht zu hoher Ordnung die 
obengenannten GréBen fiir auflésbare Untergruppen berechnet und in iiber- 
sichtlicher Form angibt. Die dabei benutzten gruppentheoretischen Satze und 
Begriffe sind wohlbekannt und recht elementar, sie werden ohne Beweise im 
2. Abschnitt zitiert, der zunachst eine allgemeine Beschreibung des Verfahrens 
gibt. Im 3. Abschnitt wird dann diskutiert, wie die speziellen Eigenschaften 
einer dualen Rechenmaschine mit hinreichend flexiblem Befehlscode zur még- 
lichst zweckmaBigen Durchfiihrung des Verfahrens benutzt werden kénnen. Im 
4. AbscHnitt wird die Art der Ausgabe der Ergebnisse angegeben. Im 5. Abschnitt 
wird eine besonders schnell arbeitende Variante des Verfahrens fiir regulare 
Permutationsgruppen beschrieben, wahrend Abschnitt 6 iiber Anwendungs- 
méglichkeiten, bisherige und geplante Rechnungen mit den vorhandenen Pro- 
grammen berichtet. In einem Anhang wird die benutzte Maschine, eine Z 22, 
kurz charakterisiert. 

Die hier beschriebenen Programme wurden in der Zeit von April 1959 bis 
Februar 1960 am Rechenzentrum der Universitat Kiel entwickelt, dessen Mit- 
arbeitern, insbesondere Herrn Dr. B. SCHLENDER, ich fiir bereitwillig gewahrten 
Rat und Unterstiitzung beim Programmieren herzlich danken méchte. 


2. Das Rechenverfahren 

Die Elemente einer zu untersuchenden Gruppe kénnen in verschiedener Form 
gegeben sein, etwa 

a) als Worte in abstrakten Erzeugenden, zu denen definierende Relationen 
vorliegen, 

b) als Permutationen, 

c) als Matrizen. 

Von diesen Méglichkeiten erschien a) wenig geeignet, da es keinen Algorithmus 
gibt, um zu entscheiden, ob zwei formal verschiedene Worte das gleiche Gruppen- 
element darstellen ; man miiBte also spezielle Reduktionen fiir bestimmte Gruppen- 
klassen einfiihren. Von den prinzipiell gleich geeigneten Méglichkeiten b) und 
c) wurde b) gewahlt, da die Multiplikation zweier Permutationen weniger arith- 
metische Operationen benétigt und, jedenfalls bei regularen Darstellungen, die 
Anzahl der pro Gruppenelement zu speichernden Ziffern geringer ist ; Rechenzeit 
und Speicherplatz unterliegen aber bei der Verwendung einer kleineren Maschine 
ziemlichen Beschrankungen. 

Da jede endliche Gruppe einer Permutationsgruppe isomorph ist, bedingt 
diese Entscheidung keine Einschrankung beziiglich der zu untersuchenden 
Gruppen. Auch fiir die praktische Rechnung treten keine Schwierigkeiten auf, 
da fiir durch Erzeugende und definierende Relationen gegebene Gruppen im 
AnschluB an die hier beschriebenen Arbeiten bereits Programme [2] hergestellt 
sind, die mit dem oben erwahnten Verfahren der Restklassenabzaéhlung Permu- 
tationsdarstellungen der gegebenen Gruppe nach beliebig vorgebbaren Unter- 
gruppen liefern. 

Es seien also n Permutationen G,, ..., G,, vom Grad g vorgegeben. Die Unter- 
suchung der von diesen erzeugten Gruppe & = {G,, ..., G,} geschieht in mehreren 


groBeren Schritten. 
21* 
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I. Die Gruppe & wird erzeugt, d.h. es wird eine Liste aller Elemente von & 
und ihrer Ordnungen hergestellt. 

II. Die zyklischen Untergruppen von & werden aufgesucht. 

III. Die minimalen Untergruppen von & werden in Klassen konjugierter 
Untergruppen geordnet. 

IV. Der Halbverband aller auflésbaren Untergruppen von © wird bestimmt. 

Diese Schritte, denen vier gr6Bere Unterprogrammsysteme entsprechen, sollen 
nun im einzelnen beschrieben werden. 

I. Um das Erzeugnis {G,,...,G,} zu erhalten, wird zunachst U,={1} und 
dann nacheinander U;= {U,_,, G,} fiir alle7=1,..., m gebildet. Dabei erhalt man 
das Erzeugnis {Ul, G} einer Untergruppe U und eines Elementes G folgendermaBen : 

Die Elemente von U seien, mit dem Einselement 1 beginnend, sonst beliebig, 


angeordnet: 
U, = 1, is sax a. 


Ist y-+4 die kleinste positive Zahl, fiir die G’+* CU, so wird B=U+UG+UG?+... 
+UG” gebildet und der Anordnung von U entsprechend geordnet: 


¥,=U,,...,V,.=U, V4.,=U,G,...,%,=UG...., %=UGC. 


Dann werden die Elemente von ¥% der Reihe nach von links mit G multipliziert: 
GV,,GV,,.... Nach jeder Multiplikation wird gepriit, ob GV;c¥. Ist dies der 
Fall, so wird als nachstes GV;,, untersucht. Ist dagegen G*=GV,¢¥% und »,+1 
die kleinste positive Zahl, fiir die G*" +? €&, so werden zu ¥ die neuen Restklassen 
UG*, UG*?, ..., UG*” hinzugefiigt und ihre Elemente im AnschluB an die schon 
vorhandenen geordnet: 


Vor+1 as Ui G*,..., Voo+1)r a U,G*, Viva) r+4 re U; G**,..., Fixemde es U, G*", 


Danach wird die Linksmultiplikation mit G bei V;,, wieder begonnen und 
das Verfahren entsprechend fortgesetzt. Man gelangt zu einem Ende, wenn bei 
der Linksmultiplikation das letzte Element von &% erreicht wird, ohne daB Ele- 
mente GV,;¢¥ gefunden werden. 

Auf diese Weise wird jedes Element von {U, G} genau einmal erhalten: Die 
Menge & besteht aus vollstandigen Restklassen nach U, keine zwei dieser Rest- 
klassen kénnen iibereinstimmen, da vor der Bildung einer jeden Restklasse 
gepriift wird, ob ihr Reprasentant bereits in einer anderen Restklasse vorkommt. 
Ferner ist nach Konstruktion mit jedem Element VE¥% auch GVE%. Daher 
erhalt man durch Induktion nach der Lange/ der Worte U,G"U,G”... U, G” 
des Erzeugnisses {U, G}, daB ¥ alle Elemente von {U, G} enthalt. Zugleich mit 
der Liste aller Elemente von & wird die Liste ihrer Ordnungen angelegt ; dabei 
erhalt man die Ordnung einer Permutation P am einfachsten als kleinstes gemein- 
sames Vielfaches der Zyklenlingen in einer Zerlegung von P in elementfremde 
Zyklen. 

II. Fiir die weitere Untersuchung der Gruppe sind die zyklischen Unter- 
gruppen von Primzahlpotenzordnung (im folgenden mit ZUPPO abgekiirzt) ein 
wichtiges Hilfsmittel. Dies soll zundchst erlautert werden: 


Im Verlauf der Rechnung werden haufig Untergruppen gebildet, die spater 
nochmals gebraucht werden. Diese miissen in einer Form notiert werden, die 
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erstens gestattet, ohne langere Rechnung zwei Untergruppen zu vergleichen, 
zweitens mit wenig Angaben auskommt. Die Angabe aller Elemente der Unter- 
gruppen wiirde sicher die zweite Bedingung schlecht erfiillen. Wird dagegen 
nur ein aus méglichst wenigen Elementen bestehendes Erzeugendensystem notiert, 
so ist im allgemeinen die erste Bedingung nicht erfiillt, da man nicht ohne langere 
Rechnung feststellen kann, ob zwei Systeme von Elementen die gleiche Unter- 
gruppe erzeugen. Man erhialt jedoch zu jeder Untergruppe U folgendermaBen 
ein Erzeugendensystem, das U eineindeutig zugeordnet ist: 

In jeder ZUPPO wird ein erzeugendes Element ausgezeichnet. Eine Unter- 
gruppe wird durch die Menge Z(U) der in ihr enthaltenen ZUPPO erzeugt, diese 
Menge ist der Untergruppe eineindeutig zugeordnet. Die Menge E(U) der aus- 
gezeichneten Erzeugenden der in U enthaltenen ZUPPO bildet daher ein U 
eineindeutig zugeordnetes System von Erzeugenden, das den genannten Ansprii- 
chen weitgehend entspricht; denn es gilt fiir zwei Untergruppen U und & einer 
Gruppe & 

- 4. UC dann und nur dann, wenn E(U) c E(&). 

2. U=% dann und nur dann, wenn E(U) = E(%). 

3. G4UG=8 dann und nur dann, wenn G1Z(U)G=Z(%&). 

Gerade diese letzte Eigenschaft erweist sich fiir die Rechnungen als sehr niitzlich. 


Im Teil II des Verfahrens werden daher die zyklischen Untergruppen auf- 
gesucht und notiert. Einzelheiten, insbesondere iiber die Art der fiir jede Unter- 
gruppe zu speichernden GréBen, werden in Abschnitt 3 besprochen, da sie weit- 
gehend auf eine duale Rechenmaschine zugeschnitten sind. 

Als nachstes soll der Verband der Untergruppen von & untersucht werden. 
Die Schwierigkeit ist hierbei, ein praktisch brauchbares (d.h. mit médglichst 
geringem Zeitaufwand durchfiihrbares) Verfahren zu finden, alle Untergruppen 
aufzusuchen. 

So scheidet z.B. die Méglichkeit aus, alle Teilmengen durchzupriifen. Man 
wird vielmehr versuchen, von schon bekannten Untergruppen ausgehend neue 
zu finden. Hierbei kann man mit den zyklischen Untergruppen beginnen und 
jeweils das Erzeugnis {Ul, G} einer schon vorhandenen Untergruppe U mit einem 
Element GqU bilden. Doch ist auch dies ein itibermaBiger Rechenaufwand, da 
man dieselbe Untergruppe auf vielerlei Weise erzeugt und erst hinterher fest- 
stellen kann, daB sie schon vorhanden ist. Man wird also anstreben, bereits von 
einer gegebenen Kombination U, G festzustellen, ob ihr Erzeugnis gleich einer 
schon vorhandenen Untergruppe & ist. Die hierfiir notwendige Bedingung Uc ¥ 
und G€& reicht nicht hin, da {U, G} gleich einer noch nicht bekannten Unter- 
gruppe von % sein kann. Um mit der Nachpriifung Uc ¥ und G €¥ auszukommen, 
mu8 man also sicherstellen, daB bei jedem Schritt der Rechnung mit einer Unter- 
gruppe auch alle in ihr enthaltenen bekannt sind. Es ist zweckmabBig, zur Be- 
schreibung dieser Situation eine neue Bezeichnung einzufiihren. 

Definition. Die.Menge der Untergruppen von @&, deren Ordnung Produkt 
von k Primfaktoren ist, werde die k-te Schicht 2, von & genannt. 

Offenbar erfiillt man die Forderung, daB mit einer Untergruppe auch alle 
in ihr enthaltenen bekannt sein sollen, wenn man die Untergruppen ,,schicht- 
weise‘‘ konstruiert. Die Schicht 2; ist bekannt, sie besteht aus den Untergruppen 
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von Primzahlordnung, die unter den ZUPPO vorkommen. Man hat also, von 
2, ausgehend 2, ,, zu konstruieren. Sei eine Untergruppe U € 2, und ein Element 
GU gegeben. {U, G} ist genau dann noch nicht bekannt, wenn Ud B oder GGB 
fiir jede schon vorhandene Untergruppe 8€2,,,. Es braucht jedoch nicht 
{ul, G}¢ 2,4, zu sein. Ob dies der Fall ist, kann im allgemeinen erst festgestellt 
werden, wenn {U, G} gebildet wird. 

Ist jedoch Gim Normalisator Ny von Uenthalten, so ist genau dann {U, G}€ 2,4, 
wenn es eine Primzahl # gibt, fiir die G?€ U. Bei dem vorliegenden Programm- 
system werden nun nur solche Elemente G € Ny zur Erweiterung von U benutzt. 
Da dann U Normalteiler vom Index # in {U, G} wird, beschrankt man sich damit 
auf das Aufsuchen auflésbarer Untergruppen von &. Auf die nicht auflésbaren 
Untergruppen konnte hier um so eher verzichtet werden, als die ersten einfachen 
Gruppen zusammengesetzter Ordnung die Ordnungen 60, 168, 360 haben, Kapa- 
zitat und Rechengeschwindigkeit der benutzten Maschine aber die Ordnung der 
behandelbaren Gruppen auf etwa 300 beschranken. 

Das benutzte Verfahren arbeitet im einzelnen folgendermaBen: 

III. Man beginnt mit der Neuordnung der minimalen Untergruppen in Klassen 
Konjugierter. Dazu wird die erste der Untergruppen von Primzahlordnung in 
der Reihenfolge der ZUPPO herausgegriffen, sie und ihre Konjugierten in eine 
neue Liste eingetragen und die eingeordneten Gruppen in der alten Liste ge- 
strichen. Unter den verbleibenden wird dann wieder die erste Untergruppe von 
Primzahlordnung gesucht usw., bis alle Untergruppen von Primzahlordnung 
behandelt sind. Dabei erhalt man die Konjugierten einer Untergruppe U, indem 
man den Normalisator Ny bildet, © in Restklassen nach Ny zerlegt und aus 
jeder Restklasse einen Repradsentanten wahlt, mit dem man WU transformiert. 

IV. Bei der Bildung von 2, ,, aus 2; kann man die Rechnung weiter abkiirzen, 
indem man die Einteilung in Klassen konjugierter Untergruppen zu Hilfe nimmt. 
Es gilt namlich: 

1. Jede auflésbare Gruppe &,,,€ 2, enthalt mindestens eine Gruppe &,€ 2, 
mit &, <1 G,4,- 

2. Ist ©, IG, 41, G,41:G,=p und H,=X1G,X, Hyyi—X1G,41,X, so ist 
Dr I Dera und $y 41: D,=2. 

Man erhalt daher alle auflésbaren Gruppen aus 2,,,, indem man fiir die 
Reprasentanten # jeder Klasse konjugierter auflésbarer Untergruppen aus 2, 
alle Obergruppen aus 2,,, aufsucht, in denen ® normal ist und dann zu jeder 
von diesen alle ihre Konjugierten bestimmt. Man muB also zu einer gegebenen 
Untergruppe ® alle Untergruppen ®* auffinden, fiir die 


(*) 1RIR*, 2. KRR=—p 
gilt. Jeder Reprasentant X einer erzeugenden Restklasse XR der zyklischen 
Faktorgruppe ®*/® hat die Eigenschaften 

(**) 14. XER, 2. XPER, 3. X7RX=—R, 4. {X, KR} = K*. 


Ist die Ordnung von X gleich p*7, «21, (p,7)=1, so hat X* die Ordnung 7, 
und aus X?€ Rt folgt X**E R. 

Ist (A, p)=1, so hat X*" die Ordnung *, und da {X?*, X*"}={X}, hat mit 
X auch das Element X*’ die Eigenschaften (**). Man erhalt daher schon alle 
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auflésbaren Gruppen ®* mit (*), wenn man zur Erweiterung von ® nur die 
Erzeugenden der ZUPPO heranzieht. Ist eine solche Gruppe #* gefunden, so 
werden ihre Konjugierten in derselben Weise wie bei der untersten Schicht 
bestimmt. 

Sind alle Reprasentanten von Klassen konjugierter auflésbarer Untergruppen 
aus 2, nach diesem Verfahren behandelt, so liegen damit auch alle auflésbaren 
Gruppen aus 2,,,, nach Klassen Konjugierter geordnet, vor, und der Vorgang 
kann mit 2, ,, wiederholt werden, um 2, , zu erhalten. 

Das Aufsuchen der Untergruppen wird abgebrochen, wenn die Schicht 2, _, 
(s = Anzahl der Primfaktoren der Ordnung von ) vollstandig behandelt ist. 


3. Zur Programmierung 


Das oben beschriebene Verfahren benutzt eine Anzahl von Rechenvorgangen 
einfacher Art wie etwa die Multiplikation zweier Permutationen, aus denen 
kompliziertere, wie etwa das Aufsuchen des Normalisators einer Untergruppe, 
aufgebaut sind. Dementsprechend besteht das vorliegende Programmsystem aus 
(insgesamt 63) ineinander geschachtelten Unterprogrammen, die, zum Teil mittel- 
bar, von einem Hauptprogramm aufgerufen werden. Man kann diese Unter- 
programme in 5 Klassen einteilen: 

1. Arithmetische Hilfsprogramme. 

2. Grundprogramme zum Rechnen mit Permutationen. 

3. Grundprogramme zum Rechnen mit Hilfsgr6éBen. 

4. Gruppentheoretische Konstruktionen. 

5. Oberprogramme I bis IV. 

Die unter 1. fallenden Programme betreffen das Rechnen mit ganzen Zahlen, 
wie Multiplikation, Division mit Rest, Primfaktorzerlegung, Bildung des gréBten 
gemeinsamen Teilers mit dem Euklidischen Algorithmus, Bildung des kleinsten 
gemeinschaftlichen Vielfachen u.a., die im Laufe der Rechnung gebraycht werden. 

Unter 2. sind die Unterprogramme zum Lesen und Drucken, Umspeichern, 
Vergleichen, Multiplizieren, Transformieren und Potenzieren von Permutationen 
gemeint. Hier war zundchst zu entscheiden, wie die Permutationen gespeichert 
werden sollten. Wahrend man bei der Rechnung per Hand meist die Darstellung 
einer Permutation als Produkt ziffernfremder Zyklen (Zyklenschreibweise) bevor- 
zugen wird, ist es in der Maschine zweckmaBiger, die Schreibweise als Abbildung 


(' a zu wahlen, da hier keine Klammersetzung mitgespeichert werden muB. 


Wahlt man bei einer Gruppe vom Grad g als permutierte Objekte die natiirlichen 
Zahlen 1,...,g, so kann man sich auf die Speicherung der Zahlenreihe i P be- 
schranken, da 7 durch die Stellung in der Reihenfolge angegeben wird. Fiir kleine 
Grade hat jede dieser Zahlen7P nur wenige Binarstellen (Bits). Um den vor- 
handenen Speicherplatz méglichst gut auszunutzen, werden daher jeweils mehrere 
Zahlen hintereinander in dieselbe Zelle gespeichert. Die Anzahl der Zahlen pro 
Zelle und die Anzahl der Zellen pro Permutation werden vor Beginn der Rechnung 
aus dem anzugebenden Grad g und der Zellenlange von 38 Bits von einem Unter- 
programm ermittelt, das auch eine Reihe von Schiebebefehlen zum_,,Hervor- 
holen‘‘ der einzelnen Zahlen aus einer Zelle herstellt (Fig. 1). 
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Die Multiplikation zweier Permutationen ist in der Schreibweise ( P) sehr 


einfach zu programmieren: Das Bild 1 PQ steht in Q an der 7 P-ten Stelle. Um 
aber bei einer Multiplikation nicht fiir jede Zahl 7 ausrechnen zu miissen, in wel- 
chen Bits welcher Zelle der Permutation Q die Zahl 1 PQ steht (dazu wire eine 
Division mit Rest erforderlich), wird vor der Multiplikation die Permutation Q 
, ausgeschoben“, d.h. die Zahlen7Q auf g feste Zellen verteilt und spater das 
Resultat zur Speicherung wieder ,,zusammengeschoben“. 


Die Ausgabe von Permutationen kann durch Unterprogramme wahlweise in 
1 
der Form (; p 

Die unter 1. und 2. erwahnten Unterprogramme geniigen bereits, um den 
Teil I des Verfahrens zu programmieren. 

Bevor die nachsten Teile des Rechenverfahrens im einzelnen beschrieben 
werden kénnen, muB noch einiges iiber die Speicherung von Untergruppen gesagt 
— he pais werden. Es war bereits unter II auf 

A—— nN ; die Rolle der ZUPPO hierfiir hinge- 
yr) 1? | 2h | oe | we | 5 | oe | wiesen worden. Legt man eine Liste 
aller ZUPPO an, so geniigt es zur 
ao) 2 | 92 | we | we | we | we | we || >sZemo Charakterisierung einer Untergruppe 
U, die Listennummern aller in U1 ent- 


) oder in Zyklenschreibweise erfolgen. 



























































ll aw iw liwinaw lw haltenen ZUPPO zu kennen. Hier ist 
nun bei Benutzung einer Dualma- 

Fig. 1. Beispiel der Speicherung einer Pe vom echine bined sche elegante Form der 
Grad 19 Speicherung méglich. Sind z ZUPPO 


vorhanden, so wird jeder Untergruppe 
eine z-stellige Dualzahl zugeordnet, deren k-tes Bit 1 oder 0 ist, je nachdem, ob die 
k-te ZUPPO in U enthalten ist oder nicht. Wir wollen diese, U eineindeutig zuge- 


ordnete Dualzahl die Kennzahl K(U) nennen. Sie wird in = |+4 Zellen ge- 


speichert. Diese Kennzahlen dienen aber nicht nur zur platzsparenden Speiche- 
rung von Untergruppen; die in dualen Rechenmaschinen vorhandenen Befehle 
zur Intersektion und Disjunktion von Zelleninhalten (s. Anhang) erlauben es 
auch, mit ihnen in einfacher Weise zu rechnen. Es ist K(U)aK(¥)=K(UO®) 
und da genau dann Uc¥, wenn UnVY=—U, kann man mittels der Kennzahlen 
mit wenigen Befehlen feststellen, ob UC% gilt. Der Disjunktion zweier Kenn- 
zahlen entspricht zwar nicht das Erzeugnis der zugehérigen Untergruppen, doch 
kann etwa die haufig vorkommende Frage, ob {U,¥}C WW durch die Abfrage 
(K(U) v K(&)) aK() 2. K(U)v K(&) entschieden werden. 

Im allgemeinen geniigt es daher, von einer einmal gefundenen Untergruppe 
nur ihre Kennzahl zu speichern. Eine Ausnahme ist bei denjenigen Untergruppen 
gemacht, die als Reprasentant einer Klasse Konjugierter auftreten. Wie in 
Abschnitt 2 erlautert, wird unter gewissen Bedingungen aus einem Reprdsen- 
tanten #¢L, und einem weiteren Element G eine neue Gruppe {R, G} erzeugt. 
Dazu braucht man aber eine vollstandige Liste der Elemente von 8. Diese 
nach den Angaben aus K() mit Hilfe des Oberprogrammes I herzustellen, ist 
wegen der vielen Linksmultiplikationen und Vergleiche in diesem langwierig. 
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Es ist zweckmaBiger, sich zu erinnern, daB ® durch sukzessive Erweiterungen 
von Primzahlindex entstand. Kennt man die jeweils neu hinzugekommenen 
Elemente X;,so kann man & leicht durch wiederholte Bildung von ¥,_,+R;_,X; 
+ ---+,_,Xf-? bei bekanntem #; ohne Abfragen herstellen. Als X; konnten 
nach Abschnitt 2 Erzeugende von ZUPPOs gewahlt werden; deren Nummern 
(ebenfalls mehrere in einer Zelle) werden daher fiir jeden Reprasentanten ge- 
speichert. (Merkzahlen M(R) fiir den Reprasentanten.) Eine Speicherung nach 
Art der Kennzahlen ist hier nicht méglich, da es auf die Reihenfolge der X; 
ankommt. 

Die zum Rechnen mit den Kennzahlen und zum Auswerten der Merkzahlen 
notigen Unterprogramme sind unter 3. zusammengefaBt. 

Das Oberprogramm II stellt die fiir die Rechnung mit diesen HilfsgréBen 
nétigen Listen iiber die zyklischen Untergruppen von & zusammen. Von I her 
liegt dazu eine Liste A, aller Elemente von & und eine Parallelliste A, ihrer 
Ordnungen vor. Nun wird zunachst die OrdnungG von & in “rimfaktoren 
zerlegt. 

G= pt py... PS,  Pi<pe<--<f,. 
Dann werden die ZUPPO in der Reihenfolge 


Pir Pes «++, PRY, Das ones PRM, 0005 Pgs eee, PPO 
(p? Exponent einer #,-Sylow-Gruppe von ®) aufgesucht und fiir jede ZUPPO {P} 
eingetragen: 

1. ein ausgezeichnetes erzeugendes Element P in eine Liste A, der Erzeugenden 
von ZUPPO, 

2. die Kennzahl K({P}) in eine Parallelliste A, zu A, sowie, um schon be- 
handelte Elemente zu kennzeichnen, 

3. in der Liste A, der Ordnungen an der Stelle der ausgezeichneten Erzeugenden 
die Adresse derselben in Liste Ag, 

4. an die Stelle der Ordnungen der zur Ordnung von P teilerfremden Potenzen 
von P eine Null. 

Die verschiedenen ZUPPO werden gefunden, indem in der Liste A, unter 
den noch nicht abgeanderten Ordnungen der Reihe nach nach Elementen der 
Ordnungen 4,, f?, ..., £%, Do, ..-, P%, D,, ---, P™* gesucht wird und fiir jedes ge- 
fundene die Operationen 1. bis 4. ausgefiihrt werden, bevor weitergesucht wird. 
Da der Exponent der #,-Sylow-Gruppe kleiner sein kann als ihre Ordnung, wird 
die Suche nach Potenzen von #,; auch abgebrochen, wenn keine Elemente einer 
Ordnung #7‘ gefunden werden, und sogleich zur nachsten Primzahl tibergegangen. 

Nach den ZUPPO werden die restlichen zyklischen Gruppen in gleicher Weise 
gespeichert. Hierbei erhalt man die Kennzahl einer zyklischen Untergruppe 3, 
die nicht von Primzahlpotenzordnung ist, als Summe der Kennzahlen der in 8 
enthaltenen ZUPPO. Diese findet man, indem man alle Potenzen eines erzeugen- 
den Elementes von 3 in der Liste A, aufsucht und nach Angabe in der Parallel- 
liste A, feststellt, ob es sich um ein erzeugendes Element einer ZUPPO handelt. 

Die Programme III und IV benutzen gemeinsam eine Reihe von Unter- 
programmen, die gruppentheoretische Konstruktionen durchfiihren, bei denen 
sowohl mit Permutationen als mit HilfsgréBen gerechnet wird. Solche werden 
gebraucht zum 
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Erzeugen einer Untergruppe nach gegebener Merkzahl, 

Bilden von {U, X}, K ({u, X}), M({Uu, X}) zu gegebenem U und X € Ny mit XEN, 
Bestimmen von Ny und K(Ny) zu gegebenen U, 

Zerlegen von & = Ny+ Ny X,+---+NyX,, 

Bilden von K(XUX) zu gegebenen A(U) und X €G, 

Aufsuchen der maximalen Untergruppen von U. 


Bei den genannten Programmen treten keine neuen Typen von GréBen auf, 
daher will ich auf ihre im einzelnen recht komplizierte Organisation nicht naher 
eingehen. 

Um auch in den Programmen III und IV selbst méglichst weitgehenden 
Gebrauch vom Rechnen mit Kennzahlen machen zu kénnen, werden auBer einer 
Hauptliste, die die nach Klassen Konjugierter geordneten Kennzahlen aller Unter- 
gruppen enthalt, drei Nebenlisten gefiihrt, die fiir jede als Reprasentant ihrer 
Klasse ausgewahlte Untergruppe zusatzlich ihre Adresse in der Hauptliste, ihre 
Merkzahl und die Kennzahl ihres Normalisators aufbewahren. Dies beschleunigt 
vor allem das Aufsuchen der Elemente X, die mit einem gegebenen Reprasen- 
tanten REL, eine Gruppe {R, X}E€2X,,, erzeugen sollen. Wie oben erklart, 
konnte man sich fiir deren Auswahl auf die Erzeugenden der ZUPPO beschranken. 
Von den an diese gestellten Bedingungen 


1. XENg, 

2. XER, 

3. {R, X} noch nicht in 2,,, vorhanden, 

4. XPER 
lassen sich die ersten 3 durch Vergleich von Kennzahlen entscheiden, da K({X}), 
K(®), K(Ng) und K(U) fiir alle UC 2,,, vorliegen. 


4. Ausgabe der Ergebnisse 


Die von der Maschine berechneten Daten werden wahrend der Rechnung 
ausgegeben. Im Programmteil I wird die Liste A, aller Elemente von & wahl- 
weise in Zyklenschreibweise oder in der Form (; 2 gedruckt, in II die Liste A, 
der zyklischen Untergruppen, wobei fiir jede zyklische Untergruppe die A,- 
Nummern ihrer Elemente angegeben werden. Zur Beschreibung der Umordnung 
der minimalen Untergruppen in Klassen konjugierter wird in III fiir jede Unter- 
gruppe ihre neue Nummer in der Hauptliste (gekennzeichnet durch ein ’) und 
ihre alte Nummer in Liste A, ausgedruckt; fiir die vom Reprasentanten # einer 
Klasse verschiedenen Untergruppen U ferner das Element P,, fiir das P;*RP=uU 
ist. Bei den Untergruppen der weiteren Schichten werden in IV ihre Hauptlisten- 
nummern und die ihrer maximalen Untergruppen angegeben; fiir einen Repra- 
sentanten ®, € 2,,,in der Form r’ = (s’, z) auBer seiner Hauptlistennummer 7 die 
seines Normalteilers ®,, aus der er zusammen mit der ZUPPO 8, erzeugt wurde, 
fiir die restlichen Untergruppen einer Klasse wieder die transformierenden Ele- 
mente. Nach AbschluB der Rechnung werden die Hauptlistennummern der Nor- 
malisatoren von Reprasentanten aufgesucht und gedruckt. Als einfaches Beispiel 
folgt das Rechenprotokoll und der danach gezeichnete Verband der Diedergruppe 
der Ordnung 12. 
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Erzeugende 
(14 2 3) (4 5) 
(1 2) (4 5) 
Alle Elemente 
P2 (4 2 3) (4 5) 
P3 (1 3 2) 
P4 (4 5) 
Ps (1 2 3) 
P6 (4 3 2) (4 5) 
P7 (4 2)(4 5) 
Ps 23 
PQ (4 3)(4 5) 
P10 (1 2) 
P11 (2 3)(4 5) 
Fmasz 
Zyklische Untergruppen 
oO. 2 £4 
4. 2 39 
2 P8 
3 P9 
4 P10 
5 Pit 
6 P12 
7 P3 PS 
8 P2 P3 P4 PS P6 
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Verband 
0’=0 
1’=1 
2’=3 TR P8 
3’= 5 TR P9Y 
4’=2 
5’= 6 ah #7 
6’= 4 TR P9 
7=7 
8’=(0' 1) ao aes 
Y TR PS 7. 2 oC 
10° TR P9 2 © 
44’ = (0’ 7) 7’ 0’ 
12’ = (7’ 1) i a 2’ 4’ 
13;’= (7’ 2) “< 6’ 5’ 4’ 
14’ 13’ 12’ 11’ 10’ 9’ 8’ 
Normalisatoren 
N( 0’) = 14’ 
N( 1’) on 8’ 
N( 4’) = 10’ 
N( 7’) = 14’ 
N( 8’) a 8’ 
N( 11’) = 14’ 
N( 12’) = 14’ 
N( 13’) = 14’ 





Fig. 2. RechenprotoKoll und Verband fiir die Diedergruppe der Ordnung 12 
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5. Vereinfachungen fiir regulare Permutationsgruppen 


Ein erheblicher Teil der Rechenzeit des obigen Programmes wird fiir das 
Multiplizieren von Permutationen und das Aufsuchen von Permutationen in 
Listen bendtigt, wie es z.B. auftritt, wenn Kennzahlen zusammengestellt werden. 
Der Zeitbedarf beider Vorginge wachst etwa quadratisch mit dem Grad der 
Permutation. 

Fiir regulare Permutationsgruppen, d.h. trans:tive Gruppen, deren Grad gleich 
ihrer Ordnung ist, kann man jedoch die Rechenzeit wesentlich verkiirzen. In 
einer regularen Permutationsgruppe gibt es namlich zu zwei vorgegebenen per- 
mutierten Objekteni und k genau eine Permutation P mit 1P=k; d.h. jede 
Permutation P ist schon durch das Bild 1 P eindeutig bestimmt. Ordnet man 
die Liste A, nach der GréBe der Zahlen 1 P, so braucht man spater nur die erste 
Zahl der neu errechneten Permutation zu untersuchen und kann in Parallel- 
listen zu A, alle iiber dieses Element wiinschenswerten Informationen speichern. 
In allen weiteren Listen geniigt es dann, die Zahl 1 P zu speichern. Insbesondere 
wird die Rechenzeit fiir die Multiplikation zweier Permutationen vom Grad der 
Permutationen unabhiangig, da in jedem Fall nur 1 PQ nach der friiher geschil- 
derten Methode gebildet wird. 

Auch in einer ganzen Reihe von Unterprogrammen ergeben sich durch Be- 
nutzung von Parallellisten zu A, erhebliche Zeiteinsparungen, insbesondere bei 
der Bildung des Normalisators. Die Anwendung dieses Programmsystems ist 
dafiir durch den hohen Speicherbedarf der Liste A, bei der Z 22 auf Gruppen 
bis etwa zum Grad 100 beschrankt. 


6. Anwendungen 


Die beiden beschriebenen Programmsysteme umfassen je etwa 3000 Befehle, 
so daB von den etwa 7000 freien Trommelzellen der Z 22 etwa 4000 zur Speiche- 
rung von Daten fiir die spezielle Gruppe zur Verfiigung stehen. Speicherbedarf 


























Tabelle 
nonent i —_— der 
Gruppe \Ordnung| Grad —S" ao xy konj alaster Rechenzeit 
gruppen Unter- 
gruppen 
S! 24 4 16 30 11 18 min 
we 60 | 5 34 59 9 40 min 
LF (2,7) | 168 | 7 78 179 15 3,5 Std 
Dive 192 | 8 61 351 58 12,5 Std 
Gro2 192 | 8 89 469 78 22,5 Std 


und Rechenzeit fiir die Diskussion einer Gruppe in der angegebenen Weise hangen 
ziemlich kompliziert von der Struktur der Gruppe (Ordnung, Grad, Anzahl der 
ZUPPO, Anzahl der Klassen konjugierter Untergruppen u.a.) ab. In einer 
Tabelle sind einige Beispiele zusammengestellt (S* symmetrische Gruppe vom 
Grad 4, U5 alternierende Gruppe vom Grad 5). Das letzte Beispiel der Liste 
zeigt, daB wesentlich ,,gréBere‘‘ Gruppen schon aus Griinden der Rechenzeit kaum 
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noch behandelt werden kénnen. Aber auch der Speicherbedarf beschrankt die 
Ordnung der vollstandig behandelbaren Gruppen auf etwa 200—300. Mit Hilfe 
einiger verbindender Handrechnung kann man natiirlich auch zur Untersuchung 
gréBerer Gruppen diese Programme zu Hilfe nehmen, die beiden Gruppen $9. 
und @, 9. sind z.B. maximale Untergruppen einer kristallographischen Gruppe 
der Ordnung 384, deren samtliche Untergruppen gesucht werden sollten. Wird 
weniger Information iiber die von einer gegebenen Menge von Permutationen 
erzeugte Gruppe (z.B. nur deren Ordnung) verlangt, so kénnen natiirlich auch 
Gruppen héherer Ordnung untersucht werden. 

Bisher sind mit den genannten Programmen eine Reihe von Gruppen durch- 
gerechnet worden, die von speziellem Interesse waren. Ferner ist bereits be- 
gonnen, insbesondere die Variante fiir reguliare Darstellungen zur systemati- 
schen Durchrechnung der Gruppen kleiner Ordnung zu benutzen. Fiir diese 
existieren in der Literatur bisher nur eine Reihe von Aufzahlungen, die je- 
doch héchstens Erzeugendensysteme angeben und sich iiberdies zum Teil wider- 
sprechen [7], S.12. 

AuBer den hier beschriebenen sind am Rechenzentrum der Universitat Kiel 
von H. FELSCH zwei gruppentheoretische Programme entwickelt worden, von 
denen eines, wie erwahnt, aus gegebenen abstrakten Erzeugenden und definie- 
renden Relationen eine Permutationsdarstellung der Erzeugenden ermittelt und 
damit die Untersuchung von durch Erzeugende und definierende Relationen 
gegebenen Gruppen ermédglicht. 

Es ware wiinschenswert, die vorhandenen Programmsysteme weiter auszu- 
bauen, insbesondere, um Erweiterungsfragen untersuchen zu kénnen, d.h. etwa 
zu gegebenen Gruppen 9% und ¥ alle nichtisomorphen Gruppen & zu finden, die 
einen Normalteiler Jt’ mit einer Faktorgruppe G/N’ = § besitzen. Ein erster 
hierzu nétiger Schritt, der auch fiir sich Interesse hatte, ware die Bestimmung 
der Automorphismengruppe einer gegebenen Gruppe. 


7. Anhang: Eigenschaften der Z 22 


Die Z 22 der Firma Zuse KG ist eine programmgesteuerte Dualmaschine mit 
einem Magnettrommelspeicher von 8192 Zellen, deren mittlere Zugriffszeit 5 msec 
betragt, sowie 25 Magnetkernspeichern. Jede Zelle enthalt 38 Bits. Im Rechen- 
werk sind Addition (+), Intersektion (A) tind Disjunktion (v) fest verdrahtet, 
die durch folgende Tabellen definiert sind: (* Ubertrag) 


+| 01 a| 014 v | 
0|04 0/00 010 
1 | 4 0* 1/04 1)44 


Negative Zahlen werden durch das in der letzter Dualstelle um 1 erhdhte Kom- 
plement dargestellt. 

Die Befehle enthalten in den Bits 1—13 die Trommeladresse, in den Bits 14—18 
die Schnellspeicheradresse. Die Bits 19—36 sind funktionelle Bits, die bis auf 
einige Ausnahmeschaltungen unabhangig voneinander wirken, die Bits 37 und 38 
dienen zur Kennzeichnung. Jeder Befehl besteht aus einem der vier Grund- 
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befehle A (Addition), J (Intersektion), U (Umspeichern), E (Sprung) mit Zusatzen, 
die Stellen und Abfragen von Bedingungen, Indexsetzen, Verschiebungen und 
Indexzahlungen im gleichen Befehl erméglichen. Die zur Ausfiihrung eines Be- 
fehls benétigte Zeit (Wortzeit) betragt 0,3 msec. 
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Tschebyscheff-A pproximationen in kleinen Intervallen II* 
Stetigkeitssdtze fiir gebrochen rationale Approximationen 


Von 
H. MAEHLY und CH. WITZGALL ** 


Im ersten Teil dieser Arbeit [5] wurden Tschebyscheff-Approximationen durch 
Polynome im Grenzfall unendlich kleiner Intervalle untersucht. Es wurde gezeigt, 
daB die ,,Fehlerkurve‘‘ gegen ein Tschebyscheff-Polynom strebt, und es wurde 
ein Naherungswert fiir die ,, Fehleramplitude‘‘ angegeben. 

Diese Betrachtungen sollen nun auf Tschebyscheff-Approximationen durch 
rationale Funktionen ausgedehnt werden. Das ist jedoch nur mit Einschran- 
kungen méglich. So ist — es wird dies in § 5 durch Beispiele belegt — die Tscheby- 
scheff-Approximation im gebrochen rationalen Fall unter Umstanden wnstetig, 
d.h. kleine Anderungen der Funktion / und der Gewichtsfunktion v kénnen groBe 
Anderungen der Fehlerkurve 6(/, v) verursachen. 

Der vorliegende zweite Teil der Arbeit wird sich vorwiegend mit dieser Stetig- 
keits- oder Stabilitatsfrage beschaftigen. Das erfordert Aussagen iiber Inter- 
polation durch rationale Funktionen. In §%3 werden daher einige elementare 
Tatsachen aus diesem Gebiet mitgeteilt. Die Darstellung geht in einigen Punkten 
iiber das hinaus, was zu den Stetigkeitsiiberlegungen des § 4 benétigt wird. 
Es erscheint dies jedoch gerechtfertigt in Anbetracht der verstreuten und schwer 
aufzufindenden Literatur iiber rationale Interpolation. Es soll zudem eine Basis 
fiir den Vergleich von Interpolation und Oskulation geschaffen werden. Das 
Problem der Oskulation wird im geplanten dritten Teil dieser Arbeit [9] eine 
entscheidende Rolle spielen. Numerische Fragen werden nicht direkt beriihrt. 

Um Riickverweisungen zu erleichtern werden die Paragraphen in allen drei 
Teilen durchlaufend numeriert. 


§ 3. Interpolation durch rationale Funktionen 


Die Schar aller rationalen Funktionen R, die sich in der Form 





R(x) = P(x) — Pothi*t oc ob Oy 
Q (4%) Jot ¥+°°* + Im *™ 


(3.1) 
mit beliebigen reellen Koeffizienten p,;, g, darstellen lassen, bezeichnen wir mit 


R (1, m). 





* Diese Arbeit wurde erméglicht durch die Unterstiitzung des Bureau of Ships, 
US-Navy, und dessen Applied Mathematics Laboratory, David Taylor Model Basin, 
unter Contract Nonr 2406(00) mit der Universitat Princeton. 

** Jetzt Syracuse University, Syracuse 10, New York (U.S.A.) und Inst. fiir 
fiir Angew. Math. der Universitat Mainz. 
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Die Summe n:=ltm 


nennen wir den ,,Grad‘‘ der Schar #i(/, m). Wenn wir im folgenden von einer 
, Darstellung‘ von RE R(l,m) sprechen, so meinen wir stets eine Darstellung 
der Art (3.1) als Quotient zweier Polynome P und Q, fiir welche 


(3.2) GradP</ und GradQsm 


ist. Zwei Darstellungen P/Q und P’/Q’ sind nur dann gleich, wenn sie in ,, Nenner‘‘ 
und ,,Zdhler‘‘ iibereinstimmen, d.h. wenn P=P’ und Q=(Q’. 

Die Darstellung P/Q einer rationalen Funktion RE R(l, m) ist offensichtlich 
nicht eindeutig bestimmt. Man kann ja P und Q mit einer reellen Konstanten 
c=+0 multiplizieren, ohne daB sich die dargestellte Funktion andert. Dariiber 
hinaus kann zuweilen ein gemeinsamer Polynomfaktor in Nenner und Zahler 
gekiirzt oder — unter Beachtung der Gradbeschrankung (3.2) — hinzugefiigt 
werden. Zur Entscheidung, ob zwei verschiedene Darstellungen zur gleichen 
rationalen Funktion gehéren, dient das 
(3.3) Kriterium. Die Darstellungen P/Q und P/Q reprdsentieren genau dann die 
gleiche rationale Funktion, wenn 

PQ=@QP 


gilt. 

Mit Hilfe von (3.3) beweist man in bekannter Weise 

(3.4) Satz. Durch n+1=l1+m-+1 threr (regularen) Punkte ist eine rationale 
Funktion RE R(l, m) eindeutig bestimmt. 

Das Interpolationsproblem fiir rationale Funktionen ist daher wie folgt zu 
formulieren 

(3.5) Interpolationsaufgabe*. Gegeben sind paarweise verschiedene ,,Inter- 
polationsstellen“‘ z,..., 2, und die zugehérigen Ordinaten Cy,...,C,. Gesucht ist 
eine rationale Funktion RE R(l, m) mit 

R(z)=c, fir 1=0,...,a. 

(3.6) Corollar zu Satz (3.4). Es gibt hochstens eine rationale Funktion R € ®t (1, m) 
die (3.5) lést. 

Die Interpolationsaufgabe (3.5) fiihrt offensichtlich auf ein lineares Gleichungs- 
system (JAcosI [4], W. E. MILNE [6]). Die folgenden Satze sollen den Zusammen- 
hang zwischen diesem linearen Gleichungssystem (3.8) und der Interpolations- 
aufgabe klarlegen. 

(3.7) Satz. Hat die Interpolationsaufgabe (3.5) eine Lésung R, so 1st jede Dar- 
stellung P/Q von R Lésung des homogenen linearen Gleichungssystems 


(3.8) P(z;) = ¢; Q(%), $=0,...,%. 
Beweis: Trivial. 


Die Anzahl der Unbekannten in (3.8) iibertrifft die Anzahl der Gleichungen. 
Ein solches homogenes lineares Gleichungssystem hat stets nicht-triviale L6- 





* Bereits Caucuy behandelte die Interpolationsaufgabe (3.5) und gab ein Analogon 
zur Interpolationsformel von LAGRANGE an. Diese ,,Interpolationsformel von Caucuy“ 
wurde von JAcoBI, KRONECKER und Netto diskutiert (vgl. [7]). 
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sungen. Ist (P, Q) eine Lésung und ist Q=0, so verschwindet P ebenfalls, weil 
es +1 Nullstellen hat. Fiir jede nicht-triviale Lésung (P, Y) von (3.8) gilt 
somit Q=+-0. Mit anderen Worten: Jede nicht-triviale Lésung des linearen Sy- 
stems (3.8) kann als Darstellung einer rationalen Funktion gedeutet werden. 
Es gilt nun der merkwiirdige 

(3.9) Satz. Es gibt stets nicht-triviale Lésungen des homogenen linearen Glei- 
chungssystems (3.8). Alle diese Lésungen liefern Darstellungen der gleichen ratio- 
nalen Funktion. 

Beweis. Es ist nur noch die Eindeutigkeit zu zeigen. (P, Q) +0 und (P, Q) +0 
seien zwei Lésungen von (3.8). Dann gilt 


¢; P(z) Q (%) = 6 Q(%) Pla). 
Fiir c;+-0 folgt daraus 
(3.10) P(2,) Q (2) = Q(%) P(%). 
Ist c;=0, so kann man zwar nicht auf (3.10) schlieBen, dafiir 14Bt sich jedoch 
der Grad von P und P durch Abdividieren von (x —z,;) erniedrigen. In jedem 
Fall folgt somit 

PO=OP. 

Nach Kriterium (3.3) ist damit Satz (3.9) bewiesen. 

Man kann also der Interpolationsaufgabe (3.5) eindeutig eine rationale Funk- 
tion zuordnen, die sich aus jeder Lésung des linearen Gleichungssystems (3.8) 
gewinnen laBt. Wir wollen diese rationale Funktion die ,,/nterpolierende‘ zu 
(3.5) nennen. Die Interpolierende hangt von den Daten der Interpolations- 
aufgabe (3.5) ab. Dies sind neben / und m die Interpolationsstellen z; und der 
Vektor der Ordinaten 


Um die Abhangigkeit von c zu betonen, bezeichnen wir im folgenden die Inter- 
polierende mit 
_* 

Aus den Satzen (3.7) und (3.9) folgt 

(3.11) Satz. Ist die Interpolationsaufgabe (3.5) tiberhaupt lésbar, so ist die 
Interpolierende R, die Lésung. Die Darstellungen von R, entsprechen dann ein- 
eindeutig den Lésungen des linearen homogenen Gleichungssystems (3.8). 

Ist die Interpolationsaufgabe (3.5) nicht lésbar, so gibt es ,,umerreichte“‘ 
Punkte (z;,¢;) mit 


(3.12) R,(%) +¢,  (R,(z) = ce ist méglich.) 
Sind P und Q beliebige Polynome, und gilt P(z;)=c;Q(z;), so ist entweder 
P(z)/Q(%) =; oder P(z;) = Q(z;) =0. Diese Uberlegung fiihrt sofort zu folgender 
Verallgemeinerung von Satz (3.11): 

(3.13) Satz. Die Lésungen des linearen Gleichungssystems (3.8) entsprechen 


ein-eindeutig denjenigen Darstellungen der Interpolierenden R,, deren Nenner und 
Numer. Math. Bd. 2 22 
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Zihler fiir jede Abszisse z; eines unerreichten Punktes (2;,c,) gemeinsam ver- 
schwinden. 

k+1 Punkte befinden sich in ,,ausgearteter Lage‘‘, wenn sie Punkte einer 
rationalen Funktion RE R(/’, m’) mit 


U+tm <k 
sind. Ist h die Anzahl der unerreichten Punkte, so gehért die Interpolierende 
R, nach Satz (3.13) zur engeren Schar R(/—h, m—h). R, interpoliert jedoch 
die /+-m+1—h erreichten Punkte, und man hat den bekannten (vgl. W. E. 
MILNE [6], Nr. 66) 

(3.14) Satz. Ist die Interpolationsaufgabe (3.5) unlésbar, so befinden sich die 
erreichten Punkte in ausgearteter Lage. 

Eine Darstellung heiBt ,,teilerfremd‘‘, wenn Nenner und Zahler keinen Poly- 
nomfaktor gemeinsam haben. LaBt R, nur teilerfremde Darstellungen (mit der 
Gradbeschrankung (3.2)) zu, unterscheiden sich also samtliche Darstellungen von 
R, nur durch konstante Faktoren, so nennen wir R, ,,eindeutig vom Typ (1, m)*‘. 
Wir haben 

(3.15) Corollar zu Satz (3.13). Ist die Interpolierende R, eindeutig vom Typ 
(1, m), so ist sie Lésung der Interpolationsaufgabe. 

(3.16) Corollar zu Satz (3.11). Ist R, eindeutig vom Typ (l,m), so hat die 
Matrix Mi(c) des linearen Gleichungssystems (3.8) maximalen Rang 


Rang Mt(c) =n +1*. 


Der Fall Rang M(c)—="-+1 ist deshalb von Interesse, weil man dann eine 
Lésung des linearen Gleichungssystems (3.8) mit Hilfe der (m-+1)-reihigen 
Determinanten von ¥t(c) geschlossen ausdriicken kann **. 

Wir wollen jetzt zeigen, daB die Interpolierende R, stetig von c abhangt. 
Wir gehen dabei davon aus, daB die Lésungen des homogenen linearen Systems 
(3.8) stetig von c abhangen. Es kommt natiirlich darauf an, welchen Konvergenz- 
begriff wir der Schar R(/,m) zugrunde legen. Wir wahlen die gleichmafige 
Konvergenz in einem abgeschlossenen Intervall [a, b]. Der zugehérige Abstands- 
begriff ist 

|R. — Rell : = Max | R, (x) — R;(x)|. 
Im folgenden Satz werden wir unter gewissen Eindeutigkeits- und Regularitats- 
voraussetzungen nicht nur die Stetigkeit von R, nachweisen, sondern auch eine 
Schranke fiir den ,,Differenzenquotienten‘‘ 


R.— Ri mm - 
eae » wo ||c—é||:= Max |¢; —é;| 
angeben. 

* W. E. MILNE [6] Nr. 66, S. 236, gewinnt (3.16) auf anderem Wege. Seine 
Behauptung, daB im Falle von Rang M(c)<+1 eine Lésung existiert, welche nicht 
eindeutig vom Typ (/, m) ist, erscheint jedoch nicht ganz korrekt: eine Lésung muB 
nicht existieren. 

** Auch die Anwendbarkeit der Cauchyschen Interpolationsformel ist auf diesen 
Fall beschrankt. 
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(3.17) Satz. Ist R; eindeutig vom Typ (l,m) und polfrei im Intervall [a, b], 
so gibt es Konstanten H, a>0, die nur von C, den Interpolationsstellen 2, ..., 2, 
und [a, b] abhdngen, derart, daB fiir alle c mit 

lle —z|| <a 
die Aussagen gelten 
I|R. — Re|| S$ Alle —e|| 
R, 1st polfret in [a, b] 
R,(%) =c,, fir 1=0,...,n. 
Bewets. 
Ag(c),..., A;(c), Bg(c),..-, By, (c) 
seien die (+ 1)-reihigen Determinanten der Matrix M(c) des Gleichungssystems 
(3.8) mit alternierendem Vorzeichen. Die Polynome P.:=))A;(c)x‘ und 
Q,:= >, B,(c)x* sind stets eine — jedoch méglicherweise triviale — Lésung 
(3.18) (FE, Q2) 
des linearen Systems (3.8). Nun ist nach Voraussetzung die Interpolierende R; 
eindeutig vom Typ (/,m). Nach Corollar (3.16) ist also Rang M(¢)=m"-+41 und 
somit 
(Fz, Qz) 0. 

Da (P,, Q,) stetig von c abhangt, gibt es eine Konstante «,>0 so, daB fiir 
||c —¢|| Sa, auch 


(F., Qc) + 0. 
P/Q, ist dann eine Darstellung der Interpolierenden R,. Die A;(c) und B,(c) 
sind ihrerseits Polynome in den Komponenten von c. Daher gibt es (etwa nach 
dem Mittelwertsatz) eine solche Konstante M, daB fiir alle ¢ mit ||c—Z|| <a 
|A;(c) — A;(¢)| SM |e —¢|| 
|B, (c) — B,@| <M lle —a]| 
gilt. Mit der Abkiirzung 
rv: = Max(|a|, ||) 
haben wir dann 
pee |P—PIISM-[t+r4--+4] lle 
I|Q. — Qel| SM -|1+7+---+| lle —el]. 
Die Darstellung B/Q; von R; ist teilerfremd. Deshalb, und weil R; im Inter- 
vall [a, b] polfrei ist, folgt 
Min | Qz(x)| > 0. 
Wegen (3.19) gibt es daher positive Konstanten «,<a,, und L, so, daB fiir 
lle —e|| Sag 
(3.20) | 
Daher folgt aus 
Re — Rell $|] J-|| | Gell (ell 2 — Fell + all @. — @ell) 


und (3.19) die in Satz (3.17) behauptete Abschatzung. 


- mt 


<L,<o. 
Q- ’ 











22* 
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Fir ||c—¢|| <a, ist nach (3.20) R, polfrei im Intervall [a,b], Nach Satz 
(3.13) ist R, sicher dann Lésung, wenn 


y(c): = IT te 


Im allgemeinen wird das Intervall [a,b] die Interpolationsstellen z; umfassen. 
In diesem Fall folgt y(c) +0 fiir ||c —¢|| <a, aus (3.20). Ansonsten kann man 
wegen y(c) +0 und der Stetigkeit von y ein positives a,<a, finden so, daB 
y(c) +0 fiir ||c —Z|| Sag ist. Damit ist Satz (3.17) bewiesen. 


§ 4. Stetigkeit der Tschebyscheff-Approximation 
durch rationale Funktionen 


Im Intervall [—b, b] sei eine stetige Funktion / und eine positive stetige 
,,Gewtchtsfunktion’’ v gegeben. Die rationale Funktion RE R(l,m) nennen wir 
,, 1 schebyscheff-A pproximierende‘‘ oder auch kurz ,,T-Approximierende‘‘ von f im 
Intervall [—, b] mit der Gewichtsfunktion v, wenn der Ausdruck 
|"; | Sol A =+||:= R(*) —f(*) 

Mae ae 


den kleinstméglichen Wert annimmt. Die Funktion 


= d(/,0): = ——f 











| 5(x»)| = R(%9) — f (2 — f (#9) 


heiBt dann die ,,Fehlerkurve‘‘. Mit Hilfe des kleinsten Wertes x,€[—b, +4-b], 
u(%») 


fiir welchen 
R—-f 
a Serres 
gilt, erklaren wir die ,,Fehleramplitude‘‘ 
A=A(f,v):= R(%o) — f(%o) 


v (%9) 





von / mit der Gewichtsfunktion v. 

ACHIESER [1] hat den Fundamentalsatz von TSCHEBYSCHEFF auf rationale 
Approximationen iibertragen. Im Gegensatz zum Polynomfall ist hier die Anzahl 
der ,,kritischen Stellen‘‘, die zur Charakterisierung der 7-Approximierenden 
ausreicht, von der speziellen Gestalt der letzteren abhangig. Zur Formulierung 
des Satzes von TSCHEBYSCHEFF fiir rationale Funktionen bendtigen wir daher 
einen neuen Begriff. 

Ist P/Q eine teilerfremde Darstellung der rationalen Funktion RE R(l, m) 
und ist 


(4.1) 


raf fir P+0 nfs Gado, 


— oo fir P=0 
so nennen wir die nicht-negative ganze Zahl 


d =d(R) := Min (/ — I’, m — m’') 


das ,,Defizit‘‘ von R in bezug auf die Schar R(/,m). Da sich sémtliche teiler- 
fremden Darstellungen derselben rationalen Funktion nur durch einen konstanten 
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Faktor unterscheiden, ist diese Definition konsistent. d(R)-+1 ist die Dimension 
des Vektorraumes der Darstellungen (3.1) von R. Ist d(R)=0, so ist R eindeutig 
vom Typ (/, m), und umgekehrt. 

Der Satz von TSCHEBYSCHEFF lautet nun: 

(4.2) Satz. Es gibt genau eine Tschebyscheff-Approximierende RE R(l, m). 
Sie ist eindeutig dadurch bestimmt, daB die Fehlerkurve an mindestens 


l+m-+2—d(R)=n+2—d(R) 
, kritischen Stellen‘ 
— OS Hy << <q 4 SH 
den Betrag der Fehleramplitude mit wechselndem Vorzeichen erreicht. 
Mit Hilfe von Satz (4.2) kann man z.B. bestatigen, daB 


R(x) =0 
die T-Approximierende R€ #(0, 1) der Funktion 
f(x) =% 
mit der Gewichtsfunktion v(x) =1 im Intervall [—1, +1] ist. 
= 
1 
ist namlich eine teilerfremde Darstellung von R. Es folgt nach (4.1) l’=—oo, 


m'=0, also d(R)= Min(0+ 00, 1—0)=1. Die Fehlerkurve lautet 
6(x) = R(x) — f(x) =— x. 


Es gibt also zwei kritische Stellen x»=—1, x,=+1. Nach Satz (4.2) reicht das 
aus um R(x) =0 als 7-Approximierende zu erkennen. 

In Teil I dieser Arbeit [5] wurden Satze (vgl. (1.1), (1.11)) iiber die Stetigkeit 
der Tschebyscheff-Approximation durch Polynome bewiesen. Den Inhalt dieser 
Satze kénnen wir kurz so fassen 


lim 6(g, w) = (lim g, lim w). 


Es handelte sich dabei um gleichmaBige Konvergenz. Die Ubertragung dieser 
Satze auf gebrochen rationale Tschebyscheff-Approximationen bereitet nun dann 
Schwierigkeiten, wenn die Grenzfunktion /= limg weniger als n+ 2 kritische 
Stellen hat. Wir schlieBen dies aus, indem wir verlangen, daB die T-Approxi- 
mierende R von {= lim g eindeutig vom Typ (/, m) ist. Im nachsten Paragraphen 
werden wir an Hand von Beispielen zeigen, daB im Falle positiven Defizits d(R) 
die Tschebyscheff-Approximation tatsachlich unstetig sein kann. 

(4.3) Satz. Es sei im Intervall [—b, 6): 
RER(l, m) T-Approximierende von { mit Gewichtsfunktion v 
SER (1, m) T-Approximierende von g mit Gewichtsfunktion w 
A=A(f, v) die Fehleramplitude von f mit der Gewichtsfunktion v. 
Ist R eindeutig vom Typ (1, m), so gibt es Konstanten K,a>0, die nur von f, v, 


l,m und b abhangen, derart, daB fiir jede weitere stetige Funktion g und Gewichts- 
funktion w mit 


If ell +14] lu — || Se 
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und 
1 1 
ets | — Weel ll 4 
|R — S|] SK (llf— ell +14] lle — all) * 
Beweis. Wir verfahren analog zum Beweis von Satz (1.1). Wir setzen wieder 


A=0 


gilt 


voraus. Es ergibt sich dann ebenso 
(4.4) Jul + |] (We —ell + alle — mI). 


R hat nach Voraussetzung kein Defizit, besitzt also einen Satz von n+2 
kritischen Stellen x;. (Im Fall A=0 kénnen wir irgendwelche n+ 2 aufeinander- 
folgende Punkte des Intervalls [—b, 6] als kritische Stellen wahlen.) Es gilt 
infolgedessen 

i 1 
(4.5) (— 4) (R() — S(x)) = — (2+ [loll | =-|)) (lf — all + alle — #1). 
Wir beweisen nun 
(4.6) Hilfssatz. Gegeben n+ 2 Punkte 
—bS%< Hy < << % 4,5) 


und eine rationale Funktion 
R= 5 ER(L,m) mit Q(x)>oO in [—6,d)], 


welche eindeutig vom Typ (1, m) 1st. Es extstieren dann zwet Konstanten H, «>0 
derart, daB fiir jede weitere rationale Funktion S=U/[V ER (l,m) mit den Eigen- 
schaften 


sais ng in [—6, b] 
(— 1)*(R(x,) — S(x%))2—a2—a, t=0,...,.m4+1 
gilt 
|R—S||SH-a. 
U 


Bewets. Geniigt die rationale Funktion S=— den Ungleichunger 


V 
(4.8) (— 1)* (R(x) — S(x)) = (— 1)°(S) — RH) 20, 

so stimmt S mit R iiberein. Wegen Q(x;)>0 und V(x,)>0 ist namlich (4.8) 
gleichbedeutend mit 


(— 1)* (P(x) V(x) — U(x) Q(%)) 20. 
Daraus folgt jedoch (vgl. (1.8)) 
PV—UQ=0, 
* Satz (1.1) (vgl. auch [2], [3]) ist der Spezialfall von Satz (4.3) fiir m=O bis auf 
die Voraussetzung ||f—g||+|4|||v—w||S«. Diese und die weitere Voraussetzung 
es | — |} all S1 sind aber unwesentlich. Sie sind willkiirlich gewahlt um 














eine Bezugsumgebung des Funktionenpaares (/, v) abzugrenzen. 
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weil ein Polynom u-ten Grades keine +1 Vorzeichenwechsel haben kann. Der 
Fall a=0 ist hiermit erledigt. 
Ist a>0O, so betrachten wir die +2 rationalen Funktionen 


a 


»._ Xs 
Eh:= FER, m), k=0,4,....0+4, 


die durch die Eigenschaft 


(— 1)' (R(x) — Ek(x))=— a, fiir alle i +h 
oder 
(4.9) E*(x;) = R(x,)+(—1)'a fiir alle ik 


eindeutig festgelegt sind. Im Gegensatz zum Beweis von Hilfsatz (1.6), wo es 

sich um Polynome handelte, ist hier die Existenz von rationalen Funktionen E* 

mit der Eigenschaft (4.9) nicht gesichert. Da die Interpolationsaufgabe (4.9) 

jedoch fiir a=0 eine Lésung hat, namlich R, und diese Lésung eindeutig vom 

Typ (/, m) und polfrei in [— 6, b] ist, gibt es nach Satz (3.17) ein solches 

a>0, 

daB fiir alle & und a mit a<@ eine in [— d, b] polfreie Lésung E* von (4.9) existiert. 
Die rationalen Funktionen E* ,,begrenzen‘‘ die Gesamtheit der rationalen 

Funktionen S=U/V mit der Eigenschaft (4.7) in folgendem Sinn 

(4.10) Min E* (x) < S(x)< Max E¥(x) fir x€|—6),)}]. 

Wir beschranken uns darauf Min E®<S zu zeigen. Wir betrachten die Bereiche 


[— 5, x), (%1, %3), (%3, %5), «+ 


und schlieBlich als letzten 
(%,41,0] falls m gerade 


(x,, 5] falls m ungerade. 
Diese Intervalle iiberdecken [— 6, b] mit Ausnahme der Punkte %,, x3, %5,.... 
An diesen Punkten ist jedoch MinE{SS trivial. Ist nun x€ (%2;_1, %2;41), SO 
betrachten wir die Werte der Funktion 
Ei_ § 
an den Stellen . 


n={* “i reap Vi < Visa 
x fir +=2), 
Es gilt fiir «+ 27 
(— 1)* (E27 (¥;) — S(y)) = (— 1)' (R(x) - S(x)) +a20. 
Ware nun 
(— 1)? (E2? (¥2;) = S(¥2;)) ses E%) (x) — S(x)>0, 


so wiirden S und E?/ die Eigenschaft (4.8) besitzen und, weil beide in [— 6, 5] 
polfrei sind, im Widerspruch zur obigen Annahme identisch sein. Also ist 


E® (x) SS(x) fir x€ (%ej-1, Xj 41)- 
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Entsprechend zeigt man 

E®(x)<S(x) fir x¢€[—), x) 

E®(x)<S(x) fiir x€(x,4,,0] falls » gerade 

E**(x)<S(x) fiir x€(x,,6] falls m ungerade. 
Damit ist (4.10) bewiesen. Aus (4.10) folgt fiir x€[—8, b] 
Min (R(x) — Ek(x)) S R(x) — S(x)S Max (R(x) — EX(x)) 
und 
|R — S|] <Max||R — Bt. 
Nach Satz (3.17) gibt es jedoch eine solche Konstante H, daB fiir aSa@ gilt 
Max ||R —E|<H.-a. 


Da H und nur von den kritischen Stellen xy, ..., X41, der rationalen Funktion 
R und dem Intervall [—b, b] abhangen, ist Hilfsatz (4.6) bewiesen. — 


Nun gilt nach (4.5) 


(— 1)* (R(x) — S(x)) = — (2+ loll ||) (l¢ all + alle — 2) 
also nach Hilfsatz (4.6) 
|R— S|] sH(2+|I>I 


} (If — ell + alle — ~}])) 


ao 
as 
| 
Damit ist schlieBlich auch Satz (4.3) bewiesen. — 
Eine abgeschwachte Formulierung von Satz (4.3) ist 


(4.11) Satz. Ist die T-Approximierende R von f mit der Gewichtsfunktion v 
im Intervall [—, b| eindeutig vom Typ (1, m), und ist 


J 
v 














falls 





If—gll+Allv—wl|s 
2+ Ill 





lim g. =f, lim w,=v, gleichmaBig in [— b,b}, 


so folgt 
lim 4(g,, w,) = d(f,v) gleichmapig in [— 6, b}. 


Die Abschatzung (4.4) laBt sich gewinnen, ohne irgendwelche Aussagen iiber 
kritische Stellen zu beniitzen. (4.4) ist daher auch dann richtig, wenn die T- 
Approximierende R nicht eindeutig vom Typ (/, m) ist. Das bedeutet aber 


(4.12) Satz. 
Aus lim g,=f und lim w,=v gleichmafig in [—}, b] 
folgt im | A(g,, w,)| =| A(¢,»)|. 


Der Absolutbetrag der Fehleramplitude A(/, v) hangt also immer stetig von 
f und v ab. Das Vorzeichen von A(f, v) haben wir jedoch recht willkiirlich durch 


A: = 6(%p) 
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festgelegt, wo % die kleinste Abszisse ist, fiir welche die Fehlerkurve ,,anst6Bt", 
d.h. |6(x9)|=|A| gilt. Hat nun die Funktion / mehrere Satze von kritischen 
Stellen, so kann bei Annaherung an / die Abszisse x, ,,springen‘‘, und daher 
das Vorzeichen von A ,,umklappen‘‘. Diese Erscheinung tritt schon im Polynom- 
fall auf. In Satz (1.13) wurde deshalb eindeutige Bestimmtheit der kritischen 
Stellen verlangt. 

_ Diese Voraussetzung ist offensichtlich zu eng. Betrachten wir etwa die beiden 
Fehlerkurven 6, und 6, (Fig.1 und 2) fiir eine rationale 7-Approximierende 
RER(l, m) mit n=l14+m=2. R habe kein Defizit. (Solche Fehlerkurven lassen 
sich natiirlich stets realisieren, indem man mit irgendeiner rationalen Funktion 
RER(l,m) ohne Defizit {:=R—6 setzt. Nach dem Satz von TSCHEBYSCHEFF 





























t | 
6e 
AW 4 
Fig. 1, Fehlerkurve normal fiir n=2 und d=0 Fig. 2. Fehlerkurve nicht normal fiir »=2 und d=0 


(4.2) ist R dann die eindeutig bestimmte 7-Approximierende zu / mit der Gewichts- 
funktion v(x) =1.) Sowohl 6,-als auch 6, stoBen in mehr als n+2=4 Stellen 
oben oder unten an. Trotzdem besteht ein wesentlicher Unterschied zwischen 
den beiden Fehlerkurven. Im ersten Fall (Fig. 1) kénnen wir héchstens vier 
aufeinanderfolgende Punkte finden, an denen die Fehlerkurve die Fehleramplitude 
mit wechselndem Vorzeichen erreicht, wahrend es im zweiten Fall (Fig. 2) fiinf 
solche Punkte gibt. Im ersten Fall existieren unendlich viele Saétze von n+-2=4 
kritischen Stellen 


%qy << %y <4, 
aber es gilt stets 
6, (9) > 0. 


Im zweiten Fall gibt es nur zwei Sdtze von n+ 2=4 kritischen Stellen 
Ky << %Q<%y, XQ << %_y << Kz, 
(%)=%,, usw.) aber wir haben 
62(%»9) > 0, 6o(%) <0. 
Im Gegensatz zum zweiten Fall gilt 
A= 6; (%9) 


fiir jeden Satz von kritischen Stellen. Hier wird man also keine Unstetigkeit 
von A erwarten. 
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Allgemein nennen wir die Fehlerkurve 4(/,v) ,,normal‘‘, wenn entweder 
6(f, v) =0 ist, oder fiir jeden Satz von kritischen Stellen —bS%<%,<:-- 
<%q41-ab 

A(f, v) —_ d(f, v) (Xo) 
gilt. Der Beweis des folgenden Satzes (4.13) unterscheidet sich nun so wenig 
vom Beweis des Satzes (1.13) in Teil I, daB er dem Leser iiberlassen werden kann. 

(4.13) Satz. Ist die T-Approximierende R von f mit der Gewichtsfunktion v 
im Intervall [—b,b] eindeutig vom Typ (l,m), ist ferner die Fehlerkurve 6(f, v) 
normal und ist schlieBlich 


lim g, =f, lim w,=v  gleichmafig in [— b,b], 


so folgt 
lim A(8;, w,) = A(t»). 


Zusammenfassend stellen wir fest: Ist die T-Approximierende R der Grenz- 
funktion / eindeutig vom Typ (/, m), dann gelten fiir die Tschebyscheff-Approxi- 
mation mit Funktionen aus # (/, m) dieselben Stetigkeitssatze wie fiir die Approxi- 
mation durch Polynome. Fallen fiir die 7-Approximierende R der Grenzfunktion 
einige kritische Stellen aus, so sind die Stetigkeitssatze falsch (Beispiele im folgen- 
den §5). Der Fall, daB R zwar nicht eindeutig vom Typ (/, m) ist, jedoch die 
normale Anzahl von kritischen Stellen aufweist, wurde nicht untersucht. 


§ 5. Beispiele 
Wir werden zwei Beispiele zur Tschebyscheff-Approximation durch gebrochen 
rationale Funktionen diskutieren. Das erste Beispiel behandelt die Schar der 
Geraden 
x+a, «areell, 


das zweite Beispiel die Schar der Parabeln 
(x = > y « reell, 


beide im Intervall [—1, +1] mit der Gewichtsfunktion v(x)=1. Im ersten 
Fall approximieren wir durch rationale Funktionen aus (0, 1) 
ro 

“ae Jorn * ’ 

im zweiten Fall durch Funktionen aus (1, 4) 
(x —£)~ Pate 
2 Jot h * 

Der erste Fall ist recht einfach, zeigt aber bereits das in § 4 erwahnte Unstetig- 
keitsverhalten. Auch im zweiten Fall ist die Tschebyscheff-Approximation un- 
stetig. Es iiberrascht, daB man hier die 7-Approximierende und die kritischen 
Stellen noch explizit angeben kann. 


Anstatt das Gleichungssystem 
R(x;) — f(x) =(—1)'Av(x), +=0,...,.m+14 
R'(x,) — f'(%) =(—AiAv(x),  §=4,.0.,0 
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fiir die 7-Approximierende R, die Fehleramplitude 4 und die kritischen Stellen x; 
direkt zu lésen (vgl. auch WENZL [8]), haben wir einen Umweg gewahlt. In 
Beispiel 1 schneiden wir die Hyperbel 1/u mit einer beliebigen Geraden cy+c, 
(Fig. 3) und bestimmen das 4 
Maximum A der Differenz 


A(u):=c+oqu——. 


Sind a und b die (positiven) 
Stellen, an denen die Diffe- 
renz A den Wert —A an- 
nimmt, so ist nach dem Satz 
von TSCHEBYSCHEFF (4.2) 
1/u die T-Approximierende 
von ¢Cy+c,u im Intervall 
[a,b]. Durch eine Ahnlich- 
keitstransformation gelangen 
wir dann zur urspriinglichen Fig. 3. Zu Beispiel 1 
Aufgabe zuriick. InBeispiel2 4 
betrachten wir entsprechend 
alle Parabeln cy+c,u+ c.u?, | 
fiir welche die beiden extre- l 
malen Differenzen A(u,) und 
A(u,) (Fig. 4) den gleichen 
Betrag haben. Man kann 
die Koeffizienten dieser Pa- 
rabeln rational durch uw, und 
uw, ausdriicken. Die Be- 
stimmung der beiden weite- 
ren Stellen a und d, an denen 
der Betrag der Differenz A 
gleich |A(u,)| wird, fiihrt 
auf eine Gleichung dritten 
Grades, die jedoch eine 
Doppelwurzel hat. Verfasser 
hoffen auf diese Weise bei 
der Durchrechnung der Bei- 
spiele Arbeit gespart zu ha- 
ben. Die Resultate lauten: 


Betspiel 1. x+a, xE€[—1, 1], Gewichtsfunktion v=1. 
T-Approximierende: (Als Vorzeichen von \/1 +a? ist das Vorzeichen von « zu 
nehmen!) 


























Fig. 4. Zu Beispie 2 


fiir «+0 
R, (x) = y1+ ws 
0 fir «=0 
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kritische Stellen: 


fir a+0: *%=-—1 fir a=0: %=-—1 
% =Vita—a %=+1. 
%,=+1. 

Fehleramplitude: 


4 ae fiir «+0 


—1 fir «=0. 


4G5 
10,4 

















Fig. 5. Fehlerkurven in Beispiel 2 fiir a=0 und «=0,225 


Als Grenzwert der Fehlerkurve 6,, wenn « durch positive Werte hindurch gegen 
Null strebt, erhalt man 


lim d.(2) =| 4 fir x=1 
a—>+0 


—«x fir —1s*<1 


wahrend 
69(x) =— x fiiralle x€[— 1,1] 


gilt. Fiir «=0 ist daher die Tschebyscheff-Approximation unstetig. 
Beispiel 2: (x — -), x€[—1, +1], Gewichtsfunktion v = 1. 


T-Approximierende: 





a ——; —. @ 
1 a E Vi+o®—1— (V1+08+ als be clex 
R,={2 ° 2 ito®+x 
5 fir «=0 
kritische Stellen: 
fir «+0: %=—1 fiir a=0: %=—1 
% = 4(— Yita%—1+2) %, = 
%_=4}(— Vita? +1+a) %g=+14. 
%j3=+1. 
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Fehleramplitude: 


1 W422 & 2 

i 7 Vite : fiir «+0 
<= 
-+ fir «a=0 


Tschebyscheff-Approximation unstetig fiir «—0 (Fig. 5). — 
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Some inequalities involving the euclidean 
condition of a matrix 
By 
F. L. BAUER* and A. S. HOUSEHOLDER** 


1. Let A be a nonsingular matrix and x any vector. Let ||x|| represent the 
euclidean norm of x, 
|| || = max ||4 %|| = Omax (A) 
the euclidean bound norm of A, and 
%(A) =|[A || -||4~|| = max (A)/m'n (A), 


where O,,x(A) and o,,(A) are the least and greatest of the singular values 
o;(A)=A}(A" A). The scalar x(A) is called the euclidean condition of A (see 
[7], [2], [3] and a forthcoming paper, for a discussion of the condition based 
upon more general norms). 


2. Given any two unit vectors é and y, let 
(1) &4%y = exp(i ®) cosy, 0S gysa/2, 
and consider the matrices 
(2) A=(Aé,An), M=A"A, 
bn ( exp (i D) cos ’). 
0 exp(i®) sing 
If ¢ is a unit vector in the plane of € and and orthogonal to &, and if 


M = (é,0)" A* A(é,C), 
then 
M=Q"MQ, 


whence 


~ — 


«(M) Sx*(Q)x(M). 
But x(M)<x(A" A), and by direct calculation 
x®(Q) =x(Q" Q) = (1+ cos g)/(4 — cos p) = cot? g/2. 


Theorem I. Jf & and y are unit vectors, A= (A &, An), and A is nonsingular, 


then 
|€%¥n|<cosy, OSpSa/2 implies (xza4)!<2(A) cot g/2. 








* Institute for Applied Mathematics, Gutenberg University, Mainz, Germany. 
** Oak Ridge National Laboratory. Operated by Union Carbide Nuclear Company 
for the U.S.Atomic Energy Commission. 
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Corollary I. If & and n are unit vectors and (&,)" M(é, n) —M, where M is 
positive definite, then 
|E¥n|<cosp, O<gySa/2 implies x(M) < x(M) cot? g/2. 
3. If 
M - “ my ' 
Mz, Mo) 
then the squared cosine of the angle between Aé and Ay is 


2 My, Me, _ 1— det M 


M11 Myo mM; M22 


_ SHANA 

I|4 &]]-|]4 ml 

Let t, and 1 x(M) be the characteristic roots of M. Then 
det M = r2x(M), 


(thy 1g)? S (M1 + M2)/2 = t(1 a x(M))/2. 





Hence 
an ae Fas. 2. Ca 
4 &|] [14 nll) (1-+%(M))? (M)+1 
< (* x*(A) cot? g/2—1 y 
) 





2(A) cot? m/2+1 


in view of the monotonic character of (x —1)/(*+ 1) as a function of x for x=1. 
Let 
cot (y/2) = (A) cot @/2. 


The last squared quotient then becomes cos?y, and with 


E=||x||-* x, n=|Iy||"* », 
the result is 


Theorem II. For any two nonnull vectors x and y, and any nonsingular 
matrix A, 
|y* x| <||x|| ||y||cosp, OSp<-2/2 implies 
\(Ay)*Ax| <||A x||||A y||cosy, cot y/2=x(A) cot g/2. 
Corollary II, For any two nonnull vectors x and y, and a positive definite 
matrix M, 
ly" x| <|x|||lylleosp, 0<@<a/2 implies 
ly¥ Mx|?<x*Mxy"Mycos?y, cot? y/2=x(M) cot? y/2. 


The inequality obtained with g=2/2 is due to WIELANDT ([6] and also private 
communication). 


4. Let Y be an arbitrary rectangular matrix of linearly independent columns 
such that A Y is defined. The matrix 


H=I-—A"Y(Y"AA"Y)1Y"%A 


is hermitian idempotent annihilating A” Y. Methods of iterative projection for 
solving systems of linear equations [4] make use of a sequence of operators of 
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the form H, each operating upon the previous residual s to form the new residual 
Hs. Then, since H?=H, 


| s|Pyl|s||* =" Me xjr* Mr, 
where, with a mild change in notation, 
M=AA®, 
and 
r=As, x=[I—MY(Y*MY)1Y"]r=AHs. 


If y is in the subspace of MY, then Hs=0O, or x=O0. In any case, x*V—0, 
x is orthogonal to Y. Hence 
x4 Mo1x=rt Mx, 
and therefore, 
|| 7 s||2/||s||2 = |r? Mi x|2/(74 M47 x? M4 x). 
Now let 
cos(Y,r) =||Y(Y" Y)2Yr|l/||7|I, 
which is the cosine of the angle between 7 and the subspace of Y. Since x is 
orthogonal to Y, 


cos(%, 7) < [4 — cos?(Y, 7) }4. 
Moreover, since 
x(M%)=x(M), (A") =x(A), 

Theorem IT, or Corollary II implies: 

Theorem III. J p is defined as in Theorem II, and 

H=I-—A*Y(Y4AA*Y)1Y7A, s=Ar, 

then 

|Y¥(Y" Y)2Y"r|| > ||r||sing, O<p~<a/2 implies ||Hs|| <||s|| cosy. 
In terms of this estimate, the optimal bound is obtained when the space of Y 
contains 7, and in that case 


|| 7 s|[/|s|| S (4? (A) — 4]/[2? (A) + 1]. 
In the usual case, Y has only a single column y, and for this, 
Corollary III. Jf H=I—A*® yy" Aj||A" y||?, then 
Tyr] Shylllirl|sing, 0<@<a/2 implies ||Hs|| <|s|| cosy. 
5. Let 
ly"7|=[Ky[llI7llsingo, OSM S2/2. 
Elementary manipulation gives 


ee 4 ee. 
IIs I? lly Al] []A-*r]] 


and from Corollary III it follows that 


MyF A|[]A*r]] < sin % 
Iy* Ilr i] sin yp © 


Moreover, sin @o/sin y= [(% +1) + (x — x1) cos go]/2 decreases if sin gp increases. 
Hence 
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Theorem IV. For any nonnull vectors u and v, and any nonsingular matrix A, 
|v «| > ||v%|| ||u||sing, OS p<a/2 implies 


{oF Al| || A> lf a B® os x +21) + (x — x) cos ~]/2 
oe ley siny — LE) + — HD cos pl/2, 
where x=x(A). 
Corollary IV. For any nonnull vectors u and v, and any positive definite 
matrix M, 
|v u| > ||v"l| ||u||sing, OS pX<a/2 implies 


ut Muvt M-v sin? @ 
Liliana < Sin’? _ [F/,,2 i. 2 2 
uH y vH v = sin? y [(2? + 1) + (% 1) cos p]*/(4x?), 


where x*=x(M). 


For p=2/2, u=v, the result is the Kantorovich inequality [5]. 


6. The statement made without proof in [4] at the end is equivalent to 
Corollary III, but in somewhat different notation. Evidently the e and 7 used 
there are related to the parameters used here as follows: 


sin y = sechy, cos y= tanhy, 
#(A) = exp e = tan g/2 cot y/2. 
Elementary transformations can be made to verify that 
cos y = tanh(y + «). 


Note that the result generalizes inequality (17) of [4], the latter resulting from 
taking 7=0. 
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Recursive Programming* 
By 
E. W. DIJKSTRA 


The Aim 


If every subroutine has its own private fixed working spaces, this has two 
consequences. In the first place the storage allocations for all the subroutines 
together will, in general, occupy much more memory space than they ever need 
simultaneously, and the available memory space is therefore used rather un- 
economically. Furthermore—and this is a more serious objection—it is then 
impossible to call in a subroutine while one or more previous activations of the 
same subroutine have not yet come to an end, without losing the possibility of 
finishing them off properly later on. 

We intend to describe the principles of a program structure for which these 
two objections no longer hold. In the first place we sought a means of removing 
the second restriction, for this essentially restricts the admissable structure of 
the program; hence the name ‘‘Recursive Programming’’. More efficient use of 
the memory as regards the internal working spaces of subroutines is a secondary 
consequence not without significance. The solution can be applied under perfectly 
general conditions, e.g. in the structure of an object program to be delivered 
by an ALGOL 60 compiler. The fact that the proposed methods tend to be 
rather time consuming on an average present day computer, may give a hint 
in which direction future design might go. 


The Stack 


The basic concept of the method is the so-called stack. One uses a stack 
for storing a sequence of information units that increases and decreases at one 
end only, i.e. when a unit of information that is no longer of interest is removed 
from the stack, then this is always the most recently added unit still present 
in the stack. For example, one can construct a stack as follows: a number of 
successive storage locations are set aside for the stack and also an administrative 
quantity, the “‘stack pointer’, that always points to the first free place in the 
stack (i.e. the value of the stack pointer may be defined as the address of the 
first free location in the stack; if the stack is empty to start off with, the stack 
pointer is equal to the start address of the portion of the memory reserved for 
the stack). When an information unit is added to the stack the stack pointer 
indicates where this unit must be stored and when this has been done, the value 
of the stack pointer is accordingly increased. To remove one or more information 
units from the stack the value of the stack pointer is suitably decreased. 





* Report MR 33 of the Computation Department of the Mathematical Centre, 
Amsterdam. 
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If we mark off, on a time axis, the moments when a unit is added to or 
removed from the stack, by using an opening bracket for the addition of a unit 
and a closing bracket for its removal, then we obtain a correctly nested bracket 
structure, in which opening and closing brackets form pairs in the same way 
as they do in a normal algebraic expression involving brackets. This is closely 
related to the circumstance that we can use a stack for storing the intermediate 
results formed in the evaluation of an arbitrary algebraic expression by means 
of elementary algebraic operations. In this case our interest is always restricted 
to the most recent element in the stack. As the intermediate results are used 
only once, use of an element implies its removal from the stack. 

A simple example may be given in illustration; the successive stack locations 
are indicated by vp, 11, ¥2,... etc. The evaluation of 


A+(B—C)x(D/E +F) 


can be split up into: v9:=A; v,:=B; vg:=C; vy: =v, —v9; vg:=D; v3:=E; 
Vg 1=Vq/Ug; Ugi=F 3 vg: =Vgt Us; Vy!=VyXVq; Vg:=Vg+,; the required result 
is formed in vy. The reader will be aware that the operations used are of only 
two different types. If we denote the value of the stack pointer by k, they are: 


1. selecting a new (explicitly mentioned) number, say X, which process is 
described by v,:=X; k:=k+1; and 
2. performing an arithmetic operation, say OP, which consists of 
k:=k—1; V,-1:=0,_,0P y,. 


If we refer to type 1 by the name of the selected variable and to type 2 by the 
operator in question, then we can also record the program by means of the 


symbol sequence: 
A,B,C, —,D,E, /, F, +, x, +. 


In this description the v’s and, what is more important, specific values of the 
stack pointer no longer appear. Here it is unnecessary to specify the values of 
the stack pointer in the text, as the computer can keep track of its value during 
execution of the programm: the v’s have become completely anonymous again, 
just as anonymous as they originally were. 

The above is well known (see for instance [J]) and so elegant that we could 
not refrain from trying to extend this technique by consistent application of 
its principles?. Let us consider for a moment the operation for the selection of 
an argument, e.g. the third ‘‘v,:—C’’. This operation can be executed without 
further claim for memory space, as we assume that the numerical value of C 





1 Without doubt, in view of the vivid interest in the construction of compilers, 
at least some of these extensions have been envisaged by others but the author was 
unable to trace any publication that went further than [1]. The referee was so kind 
to send him a copy of the report ,,Gebrauchsanleitung fiir die ERMETH"“ of the 
Institut fiir Angewandte Mathematik der ETH, Ziirich, a description by HE1nz 
WALDBURGER of a specific program organisation in which similar techniques as 
developed by Professor H. RuTiIsHAusER in his lectures are actually incorporated. 
The author of the present paper thinks, however, that there the principle of recursive- 
ness has not been carried through to this ultimate consequences which leads to logi- 
cally unnecessary restrictions like the impossibility of nesting intermediate returns 
and the limitation of the order of the subroutine jump (cf. section F 44 of the report). 

23* 
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can already be found in the memory. If, instead of C, a compound term had 
occured in the expression, e.g. C=(P/(Q—R+SxXT)), then we would have 
used v, up to v; for the calculation of this subexpression, but the nett result 
of this piece of program would still be v.:=P/(Q—R+SxXT) or v.:=C. In 
other words, it is immaterial to the “‘surroundings” in which the value C is 
used, whether the value C can be found ready-made in the memory, or whether 
it is necessary to make temporary use of a number of the next stack locations 
for its evaluation. When a function occurs instead of C and this function is to 
be evaluated by means of a subroutine, the above provides a strong argument 
for arranging the subroutine in such a way that it operates in the first free places 
of the stack, in just the same way as a compound term written out in full. 


Stacked Reservations 


In the evaluation of an algebraic expression as described above, we stack 
numerical information: the next place in the stack only becomes occupied when 
we actually fill in an intermediate result, it becomes vacant as soon as its contents 
have been used. We now consider the case that this expression—which occurs 
in what we shall refer to as the (relative) main program—contains a function 
that is to be evaluated by means of a subroutine. As far as the main program 
is concerned, all variables that the subroutine introduces for its own internal 
use are anonymous, and they should be placed in the next free places of the 
stack. Within the subroutine, however, we can make distinction between three 
types of quantities. 


The parameters. We use the name ‘“‘parameters”’ for all the information that 
is presented to the subroutine when it is called in by the main program; function 
arguments, if any, are therefore parameters. The data grouped under the term 
“‘link”’ are also considered as parameters; the link comprises all the data necessary 
for the continuation of the main program when the subroutine has been completed. 
Should one wish to do so one can leave the first free place in the stack open, 
so that the subroutine can place its ‘‘function value” there. Thereafter the next 
places in the stack are used for the parameters. (In some respects it is convenient 
if all parameters occupy the same number of places in the stack; if, as in ALGOL 
60, a parameter may be given in the form of an arbitrarily complicated expression, 
the limited parameter space in the stack can refer to a point in the memory 
where further specification of the parameter is given. The amount of information 
for the link can also be regarded as being constant. Hence the amount of stack 
space used for the parameters is constant and known for every call.) 


The local variables. In general the subroutine itself is a piece of program 
in which explicit reference is made to a number of quantities introduced by the 
subroutine. Their values are no longer of interest as soon as the subroutine has 
been completed. In the course of the execution of the subroutine they can take 
on a number of different values in succession, and their values can be used more 
than once (this is the reason why they cannot remain anonymous in the sub- 
routine itself). We allow the amount of memory space occupied by these local 
variables to be dependent on one or more parameters, e.g. one of the input 
parameters can be the length of a local vector. We restrict ourselves to the 
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case that the amount of storage space occupied by the local variables becomes 
known as soon as the input parameters are given, and that it remains constant 
during that particular activation of the subroutine. (This restriction is in ac- 
cordance with ALGOL 60; from a logical point of view the restriction is not 
essential but it simplifies matters considerably.) We therefore know at the 
beginning of the subroutine, how much memory space the local variables will 
require this time, and can see to it that they will be stored in the stack immediately 
after the parameters. 

The ‘“‘most anonymous’”’ intermediate results. During the execution of the 
subroutine intermediate results that are anonymous even in the text of the 
subroutine also play a role: for, in general, expressions will have to be evaluated 
there too, and the subroutine will in its turn call in one or more subroutines 
itself. These “‘most anonymous” intermediate results are placed in the stack 
in the same way as those of the main program, but we must now begin at the 
first free place following the last local variable. 


Consequences 


The main program used the stack exclusively for intermediate numerical 
results that were formed and added to the stack once, and were later used once 
and removed. Until then there was no need to store the stack in a random 
access memory, for our interest was at all times restricted to the youngest 
element in the stack. In principle we could have used a small magnetic tape 
that would have to move one place forward in writting and one place backward 
in reading. The fact that random access is not necessary there, is a direct con- 
sequence of the fact that these stack places as such need not be explicitly men- 
tioned in the description of the computation. 

Inside the subroutine we store the most anonymous intermediate results in 
the ‘‘top”’ of the stack in just the same way. Every reference to a local quantity, 
however, implies that one is interested in a place that is situated deeper within 
the stack, and here one is interested in random access to the stack places, in 
other words, we must be able to give the places deeper in the stack some kind 
of address. The point in the stack from which the latter is available for a sub- 
routine, is handed over to the subroutine at the moment it is called in. We 
assume a static, i.e. constant description of the subroutine to be present in the 
memory: as a result, this description must be sensitive to the dynamic specifi- 
cation, as described above, of the reference point in the stack. The static re- 
ference (static addressing) of the local variables can only occur in terms of a 
fixed position with respect to the reference point. The value of this reference 
point is derived from the value of the stack pointer at the moment of the call, 
and will therefore generally vary from call to call. 


The Link 


We must now investigate which data are to be stored in the stack under the 
heading “‘link’’. For the sake of simplicity we regard our computer as consisting 
of two parts. We refer to the memory space required for the storing of the 
program and of the stack as “the memory’, and we will call the rest ‘‘the 
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arithmetic unit’’. The following considerations—suitably interpreted—apply to 
both a built-in and a programmed arithmetic unit. 

We regard the operation ‘‘v, :=C”’ of our example as an elementary operation 
of the arithmetic unit. As a result of this operation information in the memory 
has been modified, but other changes have occured too. Before the execution 
of this operation, the state of the arithmetic unit was such that the order 
“‘vg:=C”’ was the next to be obeyed, after its execution the state is such that 
the next order is to be obeyed. Part of the arithmetic unit therefore stores 
information specifying its state (it contains something equivalent to an order 
counter). If the text of the program does not specify explicitly the values of 
the stack pointer, and the arithmetic unit is therefore obliged to keep track of 
it, the value of the stack pointer is another aspect of the state of the arithmetic 
unit. One can say quite generally that every instruction is carried out correctly, 
provided that the arithmetic unit is initially in the appropriate state, and part 
of the execution of an order is the modification of the state of the arithmetic 
unit in such a way that the next order will, in its turn, be executed correctly. 
We now consider the case that the value of C cannot be found ready-made in 
the memory, but must be calculated by means of an function subroutine. As 
we are looking for an arrangement in which it is immaterial to the surroundings 
where C comes from as long as it arrives in the desired stack location, it is 
necessary that after completion of ‘‘v,:—C’’, the arithmetic unit is left behind 
in exactly the same state as if C had been transported from the memory by 
means of an elementary operation. This must hold regardless of the complexity 
of the subroutine used to calculate C. As the subroutine is in general a piece 
of program that takes full advantage of the flexibility of the arithmetic unit, 
its execution obviously implies a series of changes within the arithmetic unit, 
and we must therefore be willing to record sufficient data regarding the state 
of the arithmetic unit in order to be able to reconstruct it later. These are the 
data that must be stored in the stack under the name “‘link’”’ when a subroutine 
is called in. (It may be possible that some of the disturbances will be annihilated 
automatically, e.g. the filling of the stack, which is, in principle at any rate, 
emptied again. For such “‘reversible’’ disturbances the reconstruction requires 
no extra measures.) 

It is clear that some form of ‘‘return address” is part of the link data. What 
the link data should furthermore include depends on the number of different 
aspects of the state of the arithmetic unit, and whether their changes are in- 
trinsically irreversible. Filling the stack is a reversible process, except when 
the language permits—as ALGOL 60 does—that a subroutine under execution 
remains unfinished on account of some criterion, and is left once for all via an 
exit other than its normal “‘return’’. A function subroutine will generally be 
called in during the evaluation of an expression, so that some stack places are 
filled with anonymous intermediate results that will never be used on account 
of the unusual exit from the function routine. In that case one must be able 
to reconstruct up to which point the stack contents are still of interest. This 
facility requires an additional record in the link. 

Furthermore there exists a certain hierarchy in the information specifying 
the state of the arithmetic unit. If we regard the state of the arithmetic unit 
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between two orders, there is certain information (order counter) that sees to 
the transition to the next order. During the execution of an order, however, 
we can distinguish between different ‘‘sub-states”: the information specifying 
these substates is of no significance at the moment of transition from one order 
to the next, and there is therefore no need to record this information in the 
link, as long as the subroutine mechanism only operates between orders. The 
specification of the substates should be recorded in the link, however, as soon 
as one allows the subroutine mechanism to operate during the course of an 
elementary operation, i.e. when (for the execution of one of its sub-operations) 
a subroutine must be called in that may possibly make use of the full flexibility 
of the arithmetic unit. 

When a subroutine, say subroutine A, is activated, it must operate in the 
stack starting at a point that only becomes known at the moment of call. The 
parameters (including the link) and the local variables of A have a fixed position- 
ing with respect to each other, but their position as a whole is only determined 
at the call: the quantity specifying the position of the whole block, is called 
a parameter pointer. The program for subroutine A can only be executed correctly 
if the current value of the parameter pointer is at the disposal of the arithmetic 
unit, in other words, the parameter pointer can be regarded as one of the aspects 
of the state of the arithmetic unit, and it must therefore be recorded amongst 
the link data. Now consider the case that subroutine A calls in subroutine B. 
We call the value of the parameter pointer indicating the parameters and local 
variables of A and B respectively, PPA and P PB respectively. If subroutine B 
is called in by A, the value PPA is recorded as part of the link formed at this 
call. The place in the stack where this link is recorded is deduced from the value 
of the stack pointer at that moment and is, by definition, recorded in the para- 
meter pointer, which now assumes the value we called PPB. When we return 
from B to A the return address is found in the stack under control of PPB, 
as well as the new value P PA of the parameter pointer; finally, at the operation 
“‘return’’, the new value of the stack pointer, which suddenly decreases, can be 
deduced from the old value of the parameter pointer (in our case from P P B) 

In this process nothing forbids A from being identical with B. The sub- 
routine only has to appear in the memory once, but it may then have more 
than one simultaneous “‘incarnation”’ from a dynamic point of view: the “‘inner- 
most’’ activation causes the same piece of text to work in a higher part of the 
stack. Thus the subroutine has developed into a defining element that can be 
used completely recursively. 


Recursive Techniques and ALGOL 60 


Every procedure is regarded as a subroutine in the sense described above. 
For the sake of simplicity a block that is not a procedure can also be treated 
as a subroutine, be it that this subroutine is only called in at one point. 

One of the complications of ALGOL 60 is that not only local variables may 
be used in each block, but that explicit reference may be made in every block 
to variables that are local in a lexicographically enclosing block. When a sub- 
routine is called in, the link contains two parameter pointer values for this 
purpose. Firstly, the youngest parameter pointer value corresponding to the 
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block in which the call occurs (as described in the previous section), secondly, 
the value of the parameter pointer corresponding to the most recent, not yet 
completed, activation of the first block that lexicographically encloses the 
block of the subroutine called in. As already mentioned, the first parameter 
pointer value plays a vital role in the return at the end of the subroutine, the 
second is indispensable in localizing the global variables in the stack. As the 
second parameter pointer, by definition, points to a link in the stack, which 
in its turn contains a second parameter pointer value corresponding to the 
next enclosing block, the arithmetic unit can trace this “chain” and, in doing 
so, will find all parameter pointer values that may be necessary for localizing 
any global variable in which it may be interested. 

One can assign a so-called block number to each block, indicating the number 
of blocks which enclose it lexicographically: the main program therefore has a 
block number =0. If the program refers to a global variable it is obviously 
necessary to specify the block in which the global variable was declared; the 
block number serves this purpose and, under control of this block number, 
the arithmetic unit can find the parameter pointer value it now needs. Further, 
the reference to a global variable will specify which variable of this block is 
required: its position with respect to the parameter pointer value just found 
achieves this. The introduction of the block numbers also makes it possible 
that the arithmetic unit has immediate access to all the parameter pointer values 
it may need. They can be stored in order of increasing block number in a so- 
called ‘‘display” (e.g. a series of index registers, numbered 0, 1, 2, 3,... etc., as 
many as the maximum value of the block number). By tracing the second 
chain of parameter pointers one can, when necessary (amongst others at the 
return, in the call of a formal procedure, and at the beginning of the evaluation 
of a non-trivial formal parameter), bring the display up to date. This needs 
—and can in general—only be done for block numbers not exceeding the 
number of the block we are about to enter. It will not surprise the reader that 
the block number is to be regarded as specifying the state of the arithmetic 
unit; in consequence it has to be stored amongst the link data in the stack. 
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Error Analysis of Floating-point Computation“ 
By 
J. H. WILKINSON 


Abstract 


This paper consists of two main sections. In the first the bounds are derived 
for the rounding errors made in the fundamental floating-point arithmetic operations. 
In the second, these results are applied in the analysis of a number of computing 
techniques for the calculation of the eigenvalues of matrices. In each case the computed 
solution is expressed as the exact solution of a perturbed version of the original 
matrix and bounds are found for the perturbations. For one of the techniques, an a 
priori bound is derived for the errors in the eigenvalues themselves. 


Introduction 


1. There are a number of examples in the literature of error analysis of 
algebraic processes performed in fixed-point arithmetic, but very little seems 
to have been done for floating-point operation. It has been suggested that 
floating-point error analysis should be carried out by the computer itself, using 
a technique in which the number of significant figures in each number is auto- 
matically assessed at each stage in the computation [7]. While this technique 
has much to recommend it for some types of computation, for many algebraic 
processes it leads to conclusions which are unacceptably pessimistic. 

In this paper we describe a type of error analysis which provides a valuable 
assessment of the effectiveness of a computing technique. The general principle 
may be illustrated by considering a specific problem, the calculation of the 
zeros of a polynomial. Suppose we are calculating the zeros of the polynomial 
/(z), defined by the equation 


La 
a, 2, 


f(2)= 


oM: 


by some specific computing procedure. We shall attempt to demonstrate that 
the calculated zeros are the exact zeros of some polynomial g(z), defined by 


om? 


g(z) = D(a, + 6) 2, 

and to find bounds for (6,/a,). We shall regard the method of calculating the 
zeros as good, if the error analysis shows (6a,/a,) to be small. This approach 
was first used by GIvENs [4] in connexion with the eigenvalue problem. 

Note that a good method by this definition will not necessarily lead to 
calculated zeros which are in substantial agreement with the exact zeros. How- 
ever a calculated zero can be in poor agreement with the exact zero only if 
the latter is sensitive to small changes in the coefficients. If this is so, then we 
may regard the determination of that zero as inherently unstable and we will 
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experience difficulty with any method of solution which involves rounding errors. 
Indeed even the multiplication of the polynomial by a constant k, provided 
this involves rounding, will lead to a substantial modification of a sensitive zero. 

In general we may define our objective in analysing any particular technique 
as follows. For computation using ¢ significant figures we shall attempt to 
determine a number, ¢,, such that we can claim that the technique will determine 
correctly all those figures of the answer which are unaltered when we perturb 
the last (¢— 4) figures of the data. 

In many cases we shall wish to regard the initial data as exact and to de- 
termine bounds for the difference between the calculated solution and the exact 
solution. Occasionally our analysis will enable us to give a worthwhile “‘a priori” 
estimate of these bounds and we give one such example in § 4. More frequently 
it will not be possible to obtain bounds which do not, explicitly or implicitly, 
contain expressions which are as difficult to compute as the answers themselves. 
An example of such a problem is the solution of the set of equations 


Ax=b. 
Even if we can show that the calculated solution is the exact solution of 
(A +6A)x=b+ 6b 


and can obtain bounds for 6A and 466, we still cannot give bounds for the error 
in the solution without some knowledge of A“ or equivalent information. Never- 
theless it is reasonable to assume that a technique for which we can obtain 
satisfactory bounds for 6A and 60 is better than one for which we cannot. Often 
our result will provide a useful starting point for the determination of ‘‘a poste- 
riori’ error bounds. It should be emphasized that this “backward” type of 
error analysis is not the only type that can be successfully applied to floating 
point computation and, in any case, there are some problems to which it does 
not seem applicable. An advantage of the method is that, if the data is not 
exact, we must, in any case, consider the effect of initial perturbations. 

The first part of the paper is concerned with the development of the in- 
equalities satisfied by the rounding errors made in the fundamental arithmetic 
operations. The remainder is devoted to the application of these results in the 
analysis of specific computing techniques. Since a considerable number have 
been analysed it was necessary to make some selection. For economy of presen- 
tation we have chosen a few related techniques, all of which are concerned with 
the calculation of eigenvalues by determinant evaluation. These were selected 
because it was felt that they illustrated the main features of the method of 
error analysis and they include some of the most-accurate methods for non- 
Hermitian matrices. 

The first example is the calculation of the eigenvalues of a symmetric triple- 
diagonal matrix by the method of bisections. This is a very important process 
and it provides an illustration of the use of the analysis to obtain an a priori 
estimate of the errors in the solution. 

We next consider the fundamental limitation on the accuracy attainable in 
the eigenvalues of a non-Hermitian matrix, A, by calculating values of det (A — 1) 
and interpolating for its zeros. The general result is applied to specific methods 
of determinant evaluation of matrices of triple-diagonal and Hessenberg forms. 
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For the evaluation of determinants of Hessenberg form we have chosen HyMAN’s 
method because the analysis shows to be false a common criticism that is made 
of this method. Further it leads naturally to the consideration of the evaluation 
of determinants of Hessenberg matrices by Gaussian elimination, and shows 
that, for such matrices, interchanges (or “‘pivoting for size’) are unnecessary. 

We conclude with an analysis of a number of related matrices of Hessenberg 
form for which the eigenproblem is pathologically ill-conditioned. These have 
been used as examples by other workers [3], and the success or failure of specific 
computing techniques when applied to them has sometimes been cited as an 
argument in favour or against these techniques. We show that their effectiveness 
on these particular matrices is unrelated to their merits in general, but is mainly 
determined by the fact that the rounding errors are favourably correlated for 
some techniques but not for others. 


Fundamental inequalities for floating-point operation 

2. We shall be concerned in practice mainly with operation in the binary 
scale. For examples quoted in the text it will be necessary to employ decimal 
computation, but the extension of the results will usually be obvious. It is 
perhaps worth mentioning that error bounds for decimal computation are often 
weaker than for binary computation. 

The exact details of floating-point operation vary a little from one machine 
to another but the differences are not very important. The results we give are 
true for computation on DEUCE or ACE. In floating-point arithmetic each 
number x is represented by an ordered pair a, the mantissa, and b, the exponent, 
such that b 

x=2'a. (4) 


Here b is an integer, positive or negative, and a satisfies the inequalities 
#Sa<1 —}>a2>-1. (2) 
The number zero, has the representation a=0O and is always treated specially. 
Gn DEUCE a=b=}0 for computation in single precision floating-point arithmetic 
In general we shall assume that the number of digits available is such that } 
never exceeds capacity and that there are ¢ binary places in a. 
All the operations on ACE and DEUCE give the correctly rounded results 
with the mantissa normalised to lie in the range defined by (2). We denote 
the results of a floating-point operation on numbers x and y by the symbols 


fi(x+y), fl(xxy) and fl(x+y), 


where it is understood that x and y are standard floating-point numbers. The 
relations (3) to (6) below follow immediately from the definitions 


fl(x+y) =(x+y) (1+), (3) 
fl(x—y) =(x—y) (1+), (4) 
fl(xxy) =xy(ite), (5) 
f(x > y) = (4/y) (1+), (6) 


where 
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in each case. The last operation is not defined for y=0. On the right of equations 
(3) to (6) all operations have their exact mathematical meaning. We may write 
equation (3) in the form 


fi(x+y) =*(it¢e)+y(14+ 8) (7) 


and this relation may be described in the following terms. ‘‘The computed 
sum of x and y is the exact sum of two modified numbers x’ and y’ which differ 
from x and y respectively by at most one part in 2’.” Note that if we have a 
double-length accumulator then when cancellation occurs in the addition of two 
numbers the sum requires less than ¢ digits for its representation. The calculated 
sum is therefore necessarily exact and e=0. This is emphasized because it is 
natural to associate cancellation with a loss of accuracy. Equation (5) may be 
expressed in the form. ‘‘The calculated product of x and y is the exact product 
of x and y’ or x’ and y where x’ and y’ differ from x and y by less than one part 
in 2'.”” Alternatively we could say that it is the exact product of x(1+ «)4 and 
y(1+e)! if this were more convenient. Similar results hold for division. 


For a continued product we have 


{1(%X%_ XK 11. X Xy_) SH Hy Hq... XH, (1+ €) (1+ &) ... (1+ &,) (8) 
where 
|e] S27. 
We may write 
f(x, X %_X +++ XX) = (A+ EZ) x %Q... %, (9) 
where 
(4—2-4*< (14+ £)S(14+27-4*71. (10) 


We may therefore compute an extended product with a substantial number of 
factors and still obtain a result which has a small relative error. On ACE using 
single-precision arithmetic we have 


(1— 2-*)""" Ss (14 E) 5 (1+2-%)"™, (11) 


so that even if » is of the order of 1,000 we have a result which is in error by 
less than one part in 10". Similarly we may perform an extended sequence 
of operations involving multiplications and divisions and obtain a result with 
a low relative error. The reader will readily verify that the order in which 
multiplications are performed affects the answer, but the result given by (9) 
and (10) is independent of the order. 

The problems involved in extended summation are a little more subtle. For 
the addition or subtraction of two numbers we always obtain a result which has 
a low relative error. If, however, we add 7 numbers there are no useful bounds 
for 


[PE (x + X%y + +++ + %,)/(%1 + Hg t+: +4%,)]—1. (12) 


Indeed //(x,+%.+--:-+,) may be zero when (x%,+%,+---+4,) is not and 
viceversa. The value of the expression (12) may also depend critically on the 
order in which the floating point operations are performed. However it is still 
true that 


[E(x + Het e+ +4,) =X trgte4+ x, (13) 








Error Analysis of Floating-point Computation 323 


where each of the x; differs from the corresponding x, by a factor which is close 
to unity. 
For let us assume that 


PU (a, + Hq o++ + %,) = my (At ef?) + met ef?) +--+ 4, (1+ 6) =A,. (14) 


Then we have from equation (7) and (14) 


U(x, + %_ tes + 4%, + X41) =f1(A, + 2,41) 


(15) 
=A, (1+ 6) + %41(1+ 6) 
where | e| <2~*. 
This gives 
U(X + %_ +++ + X41) =z ( €) >) %;(4 1+ ef”) )+ %41(4 + €). (16) 
1 

From this we deduce 

(At eft) =(1tel Ate”) t=1...7 47) 


1 et =H1+e. 
Combining these results for a sequence of values of 7, we have 
(4—2-*)"*< (1+ ef) S(144+274)" or | ef | <(r —1) 2-' approximately 
(4 a yrs = (4+ + elf) ) < (1 + “yo (18) 
(= 2 


Notice that ¢!) has the same bound as ¢!) as we might expect, since the roles 
of x, and x, are indistinguishable. An immediate consequence of the above 
result is that the upper bound for the rounding errors is least when the smallest 
terms are added first, but it is seldom convenient to take advantage of this. 

A very important operation in matrix processes is the calculation of an 
inner product. From equations (5) and (18) we have 


rv) or |e|<(r+1—1)2-' approximately. 


fU (x V+ Xe Vet +++ +%,%,) = 2 HY; (4+ ef”) (19) 
where, with the same approximation as in (18) 
je] sr2* [| S(r+2-—1)2- (¢=2...7). (20) 


The evaluation of determinants of triple-diagonal matrices 


3. In the remainder of this paper we give the application of the above results 
to a number of standard computing problems. We consider first the calculation 
of the determinant of a triple-diagonal matrix, C, for which 


i C41 = Bi 41 C41 = Vi+10° 


The determinant #,, is obtained from the sequence 


=1 A=% 
MS, 95. c005 B) 21 
p, = %, Py_-1 =P, 7 : . \ 
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For exact computation #, is the determinant of the 7™ leading principal 
minor of C. We shall show that the calculated ~, are the exact determinants 
of the minors of a matrix, C’, with modified elements «;, 8; and y;. 

Let us assume that the calculated 5, 2,,..:, P,-1 are exact results for a 
matrix of order y—1 with modified elements «;, B;,y; (¢=1,2,...,7—1). We 
now calculate ~, from equation (21) using the calculated ,_, and ,_,, and 
employing floating-point arithmetic for all operations. We have 


fl (a, ?,-1) =, P,-1(1 + &) = A, (say) | e,| = aad (22) 
11(B,y, p,-2) =B,y,b,-2(4+ Ea) = B, (say) (23) 

where 
(1— 27)? S (4+ e) S(1+274* (24) 


calculated #, = /1(A — B) 


=: A,(1+ &) — B,(1+ €3) | €s| s2~ 
=a, p,-1(1 + &) (1 + €3) bed By, P,-2(4 + E2) (1 + &3) “ (25) 
If we write 


a, =a,(4 abe Sula ) 
B, =B,(1+ €2)4 (1+ €)4 (26) 
Yr =y, (1+ &)4 ; + &,)4 


then we have 


b, =4, ,-1— Br Vs P,-—2, (27) 


where calculated values of #,, p,_,; and ~,_» are implied. The calculated value 
of p, is obtained exactly from those of p,_,, ~,_», if «,, B, and y, are replaced 
by a,, B,, B,, no adjustments of the earlier «;, 8;, y; being necessary. It follows 
immediately that the calculated #, are the determinants of the leading minors 
of a matrix C’, with elements satisfying the relations 


Oy = O (1+ £) 
B, =B,(1+ Ea) (28) 
Yr = a(1+ Ey) 

where (1 — 2~')? 1+ 2,5 (14+ 278; (4—2-#<14+ 2,5 (14278. 


Error bounds for the computed eigenvalues 
of a symmetric triple-diagonal matrix 


4. We may use the result of the last section to establish strict bounds for 
the errors in the eigenvalues of a symmetric triple-diagonal matrix calculated 
by a bisection technique described below. We now have £;=y; and we may 
assume that no £; is zero since otherwise we may deal with a number of triple- 
diagonal forms of lower order. For convenience we shall assume that the a; 
and f; are converted to floating-point representation from an original fixed 
point form satisfying 


la,/S3, || $2. (29) 
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Consequently every eigenvalue, A, satisfies the inequality 
lal <2, 
since for any matrix A, the eigenvalues satisfy the relation |A| < max )’ |a; |. 
; 4s 
7 


This situation occurs in the Givens and Householder programmes on DEUCE 
and ACE [9], [12] in which a general symmetric matrix is first reduced to triple- 
diagonal form using fixed-point arithmetic. 

The method of calculating the eigenvalues is based on the fact that if for 
any trial value of A we calculate the sequence p, ?,,..., 2, for the matrix 
(C —AI), then the number of agreements in sign between consecutive members 
of the sequence is equal to the number of eigenvalues of C which are greater 
than A. This number is denoted by s(A). The result is used to locate the r-th 
eigenvalue of C in the following way. Suppose we have two numbers a and b 
such that 

s(a)<r s(b) =r; 

so that the 7-th eigenvalues lies in the interval (b, a). We now calculate s[}(a+d)]. 
If s[(a+)]=r the r-th eigenvalues lies in the interval $(a+ 9), a. 

If s[$(a+)]<r the r-th eigenvalues lies in the interval b, (a+). 

In either case the r-th eigenvalue is located in an interval of width }(a— 6). 
In (¢+-1) steps we can locate an eigenvalue in an interval of width 2~“+” (a —b). 
Now under conditions (29) we know that every eigenvalue lies between +1 and 
therefore in (¢+1) bisections we obtain the position of the 7-th eigenvalue to ¢ 
binary places. 

The comments we have just made apply to exact computation. We now 
consider the effect of rounding errors. The numbers //(«,— A) may be computed 
exactly since it will readily be verified that all trial values of A must satisfy the 
relation 

|a,—A| <1 
and these trial values require only ¢ binary places for their representation. 
Corresponding to any value of 4 the computed sequence of #, will be an exact 
sequence for a modified matrix (C’—AJ) where the modified elements satisfy 


(a, — A) = (a, — A) (1+ &) oo 
(B;)® = Br(1+ &2) (1+ é5) 
and &, &, 3 are bounded as in (22), (24) and (25). Note that the «, and #} are 
functions of A. 
From equations (30) 

a, — a, = (a, — A) {(1+ 4) (1+ és) maa 1} 

and (31) 
B, — B, =B,{(4 + &2)4 (4 + €5)* — 4} 


|a, —a,| S(1+2~)*@-1=4 
|B, —B,| S|B,| (4+ 2)? —4] (32) 
S3[(1+ 27)? — 1] =6. 


(30) 


and thence 
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The calculated sequence therefore gives exact information about the number 
of eigenvalues of C’ which are greater than 4. Now (C’—C) is a symmetric 
triple-diagonal matrix. Its diagonal elements are bounded by 4, and its off- 
diagonal elements by 6, and hence its eigenvalues are bounded by 6,+26y. 
By Lipskn’s theorem [6] eigenvalues A; of C’, and A; of C may therefore be 
ordered so that 


|A; — A;| S 6, + 26, (33) 


and accordingly the eigenvalues of C’ lie in small segments centred on those 
of C, as in Fig. 1. Now for any value of A not lying in the interval containing 
A,, the answer to the question “‘Are there less than 7 eigenvalues which are 
greater than A?” is the same for any matrix C’ as for C. For points inside the 
A, segment the answers may or may not be the same for a C’ as for C. It is 





+. = . - 3 ." = > - 
Ag As AYA; A a; 
Y Xx 
Fig. 1 


evident that these remarks remain true even if some of the intervals overlap 
or indeed even if there are severe root clusters. (Neither C or any C’ can have 
coincident roots, since it has been assumed that no £; and therefore no f; is 
zero.) For example in Fig. 1, there are clearly less than 3 roots which are greater 
than any value to the right of X and at least 3 roots which are greater than 
any point to the left of Y. 

If we consider our bisection process, then one of two things must happen. 
Either 

(i) The 7-th root is placed in the correct interval at all stages. In this case 
the r-th root is in the final interval and this interval is of width 2~*. 

(ii) There is a first stage, the k-th, at which a wrong decision is made. This 
can occur only when the evaluation of S(A) is performed at a point inside the 
r-th interval. For example in Fig. 2, if X, is in the segment surrounding 4,, 
then our technique may decide that A, is between B,_, and X, instead of X, 
and A,_,. Subsequent bisections will clearly give points X,,,, X,42, moving 
towards A, and we cannot make another false move unless or until a later X; 
also lies in the A, segment. From this we deduce that our technique will, at all 
subsequent stages, place the 7-th root between two values of which either one 
or both lies in the A, segment. A, is therefore placed in an interval of width 2~* 
the centre of which is at a point, A, satisfying 


|4,—A| SB x 2° +6,4+ 26, 
S234 (14+2-)*—-144[(14+ 2-98 — 14] 
se QE 4g gaett  g x q-t-8 
= 13x 27-*-8, 


(34) 


A result of this type was established by Givens [4] but the present analysis 
is much simpler. 
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For the symmetric matrix we have been able to pass from an error estimate 
in the form described in the introduction to a precise estimate of the errors 
in the answers. If we return to the unsymmetric matrix we cannot deduce 
a result of comparable power. However we still have the result that the 
calculated value of the determinant of (C—AJ) will be the exact determinant 
of (C’—AI) where the bounds on the elements of (C —C’) are as good as those 
for the symmetric case. Since we are giving an analysis of the eigenvalue problem 
in a more general context in § 8 we shall not pursue the point further here. 





x x——_—_-—} x 
-) Xyo1  Xev2 Xe Ar Ak-/ 
Fig. 2 


Determinants of Hessenberg Matrices 


5. We now apply our error analysis to a more difficult problem, the calculation 
of the eigenvalues of a lower Hessenberg matrix, that is a matrix A such that 


a4,;=0 721+2. (35) 
The method we shall investigate consists essentially of evaluating the determinant 
of (A — AJ) for a sequence of values of A in order to locate a zero of this function. 


We therefore consider first the evaluation of the determinant of a lower Hessen- 
berg matrix. We shall assume that 


G; 541-0 $= { to (n — 1) 
and shall write }; for the element a; ;,,. If any 6;=0 we can consider a number 
of Hessenberg matrices of lower order. 


The method of determinant evaluation we shall investigate is due to Hy- 
MAN [5]. We put x,=1 and calculate x, to x, from the relations 


Ay, % + Ayg Xp + +++ + (a,,— A) x, +b, %,4,5=0 r=1ton—1. (36) 


The value of the determinant is then given by 


det (4 — Al) = {(E4.,x,) + [ann — A) x4} (T+, (—1)"-* 37) 


If we wish to determine the zeros of det (A — A/) the constant factor ( IT b,) (—1)"-? 
r=1 


is irrelevant and we shall ignore it. We shall show that the calculated determinant 
is the exact determinant of a modified matrix (A’ — AJ) and shall obtain bounds 
for (A’ — A). 

Let us assume that the values of x, to x; are correct values corresponding 
to modified elements in rows 1 to (ry —1). We shall now derive the modifications 
which must be made in the elements of the 7-th row so that the calculated x,,, 
is an exact value for these modified elements. 

We have first 

{l(a,,— A) = (a,, — A) (1+ €,) (38) 
|e,| S27*. 

Numer. Math. Bd, 2 24 
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For the inner product we have from equation (19), (20) and (38) 
S,= fl [ayy %1 + 4,9 %_ +++ + (a,,— A) %,] (39) 
-_ 4,1 (1+ m) y + a,2(41+ m2) Xe + ate (Q.— A) (4+ é,) (1+ 7,) x 


where 
(4—2-"si1+tms(it2-y (40) 
4—2-)rt?-*§<ity,s(1+2‘)t?* s=2...7 
%.4, is now calculated from the relation 
YH = ad k S,/6,) = — S,/b,(4+¢) (41) 
where |¢| <2~‘. 
The computed x,,, is therefore derived exactly from the computed My, Xgy 005%, 


and the modified coefficients a,, to a,, and b; where 

a,,=4,,(it+tn) t=1tor—1 

(a,, — A) = (a,, — 4) (1+ €,) (1+-,) (42) 
The function of which we are finding the zeros is S,,. Unless we assume some- 
thing about the elements of A we cannot make any assumptions about the values 
of A we should use. Let us assume that the a,; satisfy the relation |a;,;| <1; 
then we can restrict A to the range | A| Sm. 


We now see that even for matrices of quite large orders the a;; are very 
close to the a;;. Suppose, for example, »=255 and ¢= 48, then since we have 


(4 — 2748)255 > 4 — 256 x 2748 — 4 — 2740 


and 
(1 + 2748)255 < 44 2-40 
and | A| <255, we can be sure that 
|a;;—4,;,| $2. 

In all cases the limits are strict upper bounds. We would expect the statistical 
distribution of rounding errors and the variation in the sizes of the a;; to ensure 
that the matrix A’, determined in the manner we have described would differ 
from A by far less than the upper bound we have obtained. For example the 
factor (1+) associated with a,, in equation (40) is the product of 7 terms of 
the form (1+ ¢;) where in general the ¢; are unrelated and all satisfy | e; <2. 
We would expect (174 2-4) to be more realistic bounds for (1+), giving 
in our example 
|a;; — a;;| S2-“4(all 7,7). 

In a recent paper [10] we examined the rounding errors made in the evalua- 
tion of a determinant by triangular decomposition with interchanges and showed 
that, in general, for fixed-point operation the calculated value is the exact 
determinant of A’ where the elements of (A’— A) are bounded by 2-“+»), The 
result is for a matrix of general form and therefore applies to Hessenberg matrices 
for which, however, the volume of computation is far less. We have considered 
here the alternative evaluation, in spite of its inferior accuracy, because it 
illustrates our technique and leads to interesting results. 








Error Analysis of Floating-point Computation 329 


Gaussian Elimination for Hessenberg matrices 


6. There is one aspect of the result of the last section which calls for special 
comment. We observe that the presence of small elements 0; is in no way harmful 
in spite of the fact that we divide by these elements. Indeed since 


b, = b, (1+ €) 


the difference between b; and b, diminishes with b,. There is no correlation 
between the accuracy of the computed determinant and the size of the b;. Further, 
if any of the }; is zero then we may replace it by a very small number ),, 2~®° 
for example on ACE, and proceed. The value obtained for the determinant 
will not suffer. 

It is quite commonly believed that small }; are harmful (see, for example 
Waite [8]) and at first sight it is surprising that this should not be so. HyMAN’s 
method is very closely related to Gaussian elimination on the transpose of the 
matrix A. Now for Gaussian elimination it is well known that the use of small 
pivots and very large multipliers usually leads to serious errors both for the 
evaluation of determinants and the solution of linear equations (see, for example, 
WILk1nson [10]). We will now demonstrate that for upper Hessenberg matrices, 
Gaussian elimination may be performed in floating point without interchanges 
for determinant evaluation, (even if it leads to very large multipliers), but not 
for the solution of linear equations. 

Gaussian elimination without interchanges of an upper Hessenberg matrix 
A® may be described as a successive pre-multiplication by matrices M,, M,,..., 


M,_,. The general form of the M; is illustrated by displaying M, for »=6. 
m™ 0 0 0 0 OF 
0100 0 0 
0010 0 0 

M,= 4 

= 10 0 ms 1 0 ~«0 (43) 
, es ys 
10 0000 4, 








The quantity m; is chosen so as to make equal to zero the (7+ 1, 7) element of 
the matrix A‘t? defined by 


A+) — M,M,_,...M,A, (44) 


The matrix A” is triangular and the product of its diagonal elements is equal 
to the determinant of A. Note that one row only of A" is modified when 
each M, is applied and that 


al) =al) k=1 to (i—1). (45) 


The multipliers m; are determined from af"), ; and the calculated a‘) from the 


relation 
— nn; = fl(ae)s > as) 


6 
= al!) s,i(1+ @)/al). - 


24* 
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(We shall use ¢ to denote any quantity, not necessarily the same each time, 
satisfying || <2-'.) The mathematical equations defining a{') and their com- 
putational equivalent are given below. In the computational equivalents all 
symbols other than the ¢ refer to actual numbers arising during the computation. 


Mathematical Computational 
al) = al!) + m, aj} af) = (1+ e) [al + m, (1+ €) af] 
af’) = af!) + m, af) af) = (1+ e) [a$!) + m, (4+ €) af] (47) 
at’) = af +m, a) al’) = (4+ e) [af!) + mg (1+ €) af] 
aS) = a) + m,_, a=", af’) = (4+ e) [a) + m,_,(1+ ) a » 
Combining these equations we have 
al’) = (1+ e) al) + m,_,(1+ e)8af)s s+ m,_1m,_o(1+ €)SaMe, +> (48) 


+ M;_1M,_9... mM (1+ €)**-§ af) roe 5-1 M;_g-.. (1+ €)**-* af), 
We conclude from equations (46) and (48) that the calculated diagonal elements 
of A™) are exactly equal to the diagonal elements which would be obtained by 
performing exact Gaussian elimination on a modified matrix A’ with elements 
in the i-th column defined by 


ai, ;= =al), i(1 + é) 
aj, =alP(1te)*Ott GF = 2... (49) 
ay; = al!) (4+ «24-9, 
The modified matrix would also give exactly the m,; obtained in practice. 
We can therefore conclude that the product of the computed diagonal elements 
of A™) is the exact determinant of A’. The bounds for (A’— A) are about twice 
as big as in HyMAN’s method and much bigger than in fixed-point operation 


using inter-changes. Note, however, that the modified matrix would not give 
the remaining elements of A™ exactly. We are therefore not entitled to make 














Table 1 

A) 6 
.12345 x 10-8 .31263 -21487 .87652 
-51065 -21053 — .31065 .43214 
-41032 — .21007 81765 

A@) c 
.12345 x 1078 -31263 -21487 .87652 

— .12932 x 10’ — .88881 x 10° —.36257 X 10” 

— .49268 — .33275 





any corresponding claim for the solution of a set of linear equations with a matrix 
of Hessenberg form using Gaussian elimination without interchanges. 

This difference in the accuracy obtained in evaluating the determinant of 
a Hessenberg matrix and solving the corresponding set of linear equations cannot 
be overemphasised. It is well illustrated in the example given in Table 1. 
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The determinant calculated from A‘ is 0.078654, which is correct to five signi- 
ficant decimals, the working accuracy. The solution of the linear equations 
calculated from A‘) x=c has however an x, component which is wrong in lla 
its figures as is shown in Table 2. 

One element only is incorrect because only on the first Table 2 
occasion was a small pivot used. In general if we do “Gichateax | Corect + 
so at the 7-th step, the first r variables will be incorrect. 











If the above example is solved by Gaussian elimination 8.6326 29260 
with interchanges then both the determinant and the 2.3395 2.3395 
solution of the linear equations are obtained correctly -67539 67539 


to working accuracy. We may carry the above con- 

clusions to their limit. If a pivotal element is zero we may replace it by 
any very small number, 107! say, and continue. The accuracy of the com- 
puted determinant will not suffer thereby. 


The calculation of the eigenvalues of Hessenberg matrices 


7. We have seen that when we attempt to compute det (A — AJ) we obtain 
the determinant of a matrix (A’—AJ) where A’ is, of course, a function of A. 
If the eigenvalues of A are A,, Ag,..., 4,, and of A’ are Aj, Aj, ..., A, then we 
have 


det (A — Al) =]] (4,—A) 


det (A’— AI) = J] (A; —A). (50) 


For computation of sufficient accuracy we know that if the A; are distinct we 


have 
(A; — A) ~ vo] (A’ — A) u](v7 4) (54) 


where v; and u; are left and right hand eigenvectors of A corresponding to 4;. 
If w; and v; are normalised so that their Euclidean norms are unity then if 


|ai;—4,;| Se, 
equation (51) gives 
|A; —A,| Sn e/|v7 u,|. (52) 


Now for HyMAN’s method we know that 


je] S(4+2-)"—-1 


; nt (53) 
—=nxX2 
for any reasonable values of n. Hence 
[Aj — aj] Sn? 27] of u,]. (54) 


The accuracy attainable is limited by the size of the inner product (v}™,) for 
normalised v; and u;. If any (v} u,) is as small as 10~*, we are almost certain 
to lose k more figures in the corresponding eigenvalue of the matrix than we 
would, for example, for a symmetric matrix, since for a symmetric matrix all 
(v7 u;) are equal to unity. We still have to examine the effect of trying to 
determine the zeros of det(A —AJ) from the computed values of (A’—AlJ). 
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Now it may well happen that some of the eigenvalues are well-conditioned 
(with the corresponding (v/ u,) near to unity) while others are very ill-conditioned. 
Since all eigenvalues are involved in equations (50) for the computed determinant 
it is pertinent to ask whether the presence of ill-conditioned eigenvalues will 
make it impossible to locate accurately even those which are well-conditioned. 
We now show that this is not so. 


Let us consider a real eigenvalue, A,, of a real matrix. Using HyMAn’s 
method we know that the computed determinant for each value of A is the 
exact determinant of some (A’(A) — AJ) and that 


|a;;—4,;| Sn2~- (55) 


for all permissible A, i and 7. Let us suppose that A, is a well-determined eigen- 
value but that some of the others are less well-determined. Suppose that for 
all possible real matrices A’ satisfying (55), the eigenvalues A; satisfy 


the 6; being different in general for each eigenvalue. By our hypothesis 4, is 
much smaller than some of the other 6;. The general case will be sufficiently 
illustrated by considering a specific example. In Fig.3 below A; is a well- 
conditioned eigenvalue and remains real for all permissible A’, A, and A, are 
a complex pair which are considerably worse-conditioned, A, and A, are an ill- 
conditioned pair of close real eigenvalues which become a complex pair for some 
permissible A’ while A, is an ill-conditioned eigenvalue which remains real. We 
have shown the domains in which each of the eigenvalues varies for permissible A’. 
Consider now the determination of the well-conditioned eigenvalue, A;, by the 
method of bisection. If we calculate the value of det(A — AJ) at points A and B 
respectively, it is clear that computed value will be negative at A and positive 
at B wherever the A; lie in their permitted regions. If then we bisect this region 
repeatedly in the obvious way, in an attempt to locate the zero, the computed 
values of the determinant will always have the correct sign unless a point of 
evaluation lies in the region surrounding A,;, when it might give the wrong sign. 
For example, we have shown the first bisection point, C, lying inside the A, 
interval. Det (A — AJ) has a negative sign at C but we may well obtain a positive 
value. In this case the next bisection point will be chosen between A and C 
instead of C and B. However, the next few bisection points will lie outside 
the 4; interval and between A and C. Computed determinants will therefore 
be negative and bisection points will move steadily towards A,; we cannot make 
another false move until a bisection point lies inside the interval surrounding A,. 
The maximum error in the answer cannot therefore be greater than (6,+% the 
final length of the bisection interval). The ill-condition of the other eigenvalues 
is of little importance. 

It may well be felt that the success is due to the fact that the bisection 
technique is unsophisticated in that it makes no quantitative use of the computed 
values. However this is not, in general, true and methods based on successive 
inverse linear interpolation for zeros or the quadratic interpolation method of 
MULLER [1/1] are no less successful. 
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We may show how this comes about by a simple hypothetical example. Let 
us assume that the matrix A of order 3 has eigenvalues 1, 2 and 3 and that for 


all matrices A’ satisfying 
—3x10° < (a;;—4,;) <3 x10°° (57) 
the corresponding eigenvalues Aj, Az, Ag satisfy 
— 108 < (A, — 1) <10°8 
— 10-6 < (A, — 2) < 10°8 (58) 
— 1072 < (As — 3) < 10°, 


so that A, is very well conditioned, A, much less so, and A, very ill-conditioned. 
Each time we evaluate det(A —AJI) for a permissible value of A we obtain 
(A —A;) (A— As) (A— Ag) instead of (A—1)(A—2)(A—3), where in general the A; 











Fig. 3 


will be different each time. We show that the determination of A, by repeated 
linear interpolation does not suffer. We calculate a sequence of approximations 
x, to the eigenvalue 1, in which x,,, is determined from x, and x,_, and the 
calculated values of /, and /,_, of det(A —AI), by the relation 


4 Up Baa — fra eM — bya) = 9 + Sema (59) 
(fr —fr—a) 
In Table 3 we compare the progress of computations starting from x,=1.2 and 
%,=1.1, for true evaluation of (x —1)(x—2)(x—3) and for (x — Aj) (x — A) 
(x — 3), where in order to simulate the effect of rounding errors we take 


A=1+10°8 ”A=2+4+20% A=3+107 


and make a random choice of the signs for each evaluation. The choice of signs 
is given against each evaluation. It is evident that both the speed of convergence 
and the accuracy finally attained are scarcely affected, although all evaluations 
are subject to errors of the order of 1 part in 200, and these errors are in random 
directions. At the beginning of the computation this obviously does not matter 
very much. When we are approaching the limit and might expect to need more 
accurate values, we find that, corresponding to the fact that the computed /, 
have errors of 1 part in 200, it is the computed values of (x,,,—x,) which have 





J. H. WILKINSON: 


334 


much the same relative error and at this stage this amounts to a very small 
relative error in the calculated x,,,. 

Much the same effect is observed with the generalised Muller technique of 
quadratic interpolation and the accuracy attainable with this technique is 




















Table 3 

x | (#1) (x2) (#—3) x signs | (x—Ai) (x—A,) (2-45) 
1.2 + .288 1.2 ++— + .2864 
1.1 +.171 1.1 —++ +.1719 
0.954 — .09845 0.950 —+— —.1071 
1.00734 +.014519 1.00758 —+— + .01491 3 
1.00048 5 + .00096 929 1.000542 4 +.0010885 
0.99999 462 — .00001 07601 0.99998 785 wee — .00002 4422 
1.00000 000 1.00000001 


astonishingly high. As examples we give the computed eigenvalues of a complex 
triple diagonal matrix of order 16 and of a complex Hessenberg matrix of order 
10. In both cases the eigenvalues were found as zeros of the A-determinants using 


Table 4. Zeros of triple-diagonal matrix 

















MULLER’S methods. The 
determinants were evaluat- 


as Bi Eigenvalues ed by the method of §3 
and HyMAn’s method re- 
3+2i +2.06853152 | —2.054430457 spectively, using  single- 
Sata] Sake] S4SURSE [HERRERA preion “touting binary 
—4t t ° oD. 4 : ° 

24311 3+ 245640400 | +0.631936861i arithmetic on DEUCE. For 
—eie 3—2i 2.27740066 | +1.44826850i convenience of presentation 
1+2i] 2-2: 0.812811959 | +1.33551135% | wechose matrices in which 
5+21 aii Wh oo PR in ol the elements were small 

—2+% 1+32] —2.724803 0.6570645467_: : 
hi 242% Peg bon - aaa integers but this does not 
—~1—4i] 343i] +3.28048252 | +3.275661631 Contribute to the accuracy 
2+i | —14+5i] +1.19252750 | —5.44399752i of the answers. In the 
1—Si}] 4+3%¢] +3.55339888 | +1.264656317 triple diagonal matrix all 
ape oar “ls a Fe seion y; were taken equal to 1.0, 

t wo &y ° 4 . ° 

~ $435 5—i | —5.65716067 | +1.63200082i1 Since the eigenvalues de- 
1—5i]—3—4i] 5.77408994 | +2.839335911 | pend only upon the pro- 


ducts $;y;. The eigen- 


values had a maximum error of 2 in the last figure retained on DEUCE, 


that is, the 30th binary figure of the mantissa. 


The analysis by this technique is so simple and automatic that it is tempting 
to think that, using floating arithmetic, one will always obtain the exact solution 
of a problem which is close to the original. This is of course not true. It is 


instructive to attempt to prove that the vector x, %,.. 


., x, calculated by 


Hyman’s process for a value of A which differs by 2~' from a true eigenvalue, 
is an exact eigenvector of a neighbouring matrix A’. No such result is in fact 
true and indeed the calculated vector may be almost orthogonal to the true 
eigenvector, although the exact solution with an exact eigenvalue, is of course 
an eigenvector. 
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The Eigenvalues of an ill-conditioned matrix 


8. The two matrices treated in the last section 
were both quite well-conditioned. The results 
confirm that when none of the (v/ u,) is small, 
eigenvalues are obtained very accurately. This 
has certainly been our experience. The method 
of HyMAN has proved to be among the most 
accurate of those we have tried, though we have 
several methods, of much the same accuracy, 
which are much faster. We now consider the 
calculation of the zeros of a family of ill- 
conditioned matrices, A,,, which has been used 
as a test matrix by others [3]. The matrices 
are defined by 








1 1 "7 
122 

A, =| 1 2 3 3 
1234... (n—1) (n—1) 
11234 m—1) on 


and we consider A,,. It will readily be verified 
that the det A, =1 for all m. The value of the 
determinant is, however, extremely sensitive to 
changes in the elements near the bottom left 
hand corner. If we denote by A(z,7) the matrix 
A,» with its (7,7) element increased by 6 then 
we may prove quite trivially that, 


det A(12,1)=1—6x11! -=1—4.00x1076 
det A(1,2) =1—6 

det A(12,2)=1+6x11! +1+ 4.00x1076 
det A(41,1) =1+612!/41 = 1+ 4.35 x1076 


det A(1,4) =1+20. 


In general the sensitivity decreases as we move 
towards the top right hand corner. Now the 
determinant of a matrix is the product of its 
eigenvalues. Therefore if we consider matrices 
Ai differing from A,. by quantities of order 
10~®, it must be true for some Aj, that at least 
some of the eigenvalues differ substantially from 
those of A,,. To emphasize this we give in 
Table 6 accurate eigenvalues of A and of A(12,1) 
with 6=10°*. It will be seen that the smaller 
eigenvalues are changed beyond recognition, 


Table 5. Zeros of Hessenberg matrix 
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while the larger ones are scarcely affected. Corresponding to this, we find that 
the values of (v/ u;) for normalised v; and u; are very small for the smaller 
eigenvalues, but of the order of unity for the larger ones. 


Table 6 





Accurate eigenvalues of 


Eigenvalues computed by Hyman’s method 
(30 bits mantissa) 

















Ais A(12,1) for 6=10-° ie | Bis 
.03102 80606 .02563 19733 .03102 80549 1st two eigenvalues 
.04950 74292 .06665 42393 .04950 74292 unrecognisable 
.08122 76592 +.01167 405232 .08122 76573 .08291... 
-14364 65198 -14677 19886 -143646520 x. Se 
.28474 97206 .28443 79091 .284749720 oy ee 
.64350 53190 .6435144692 -64350 5319 643511... 

1.5539887091 1.5539885977 1.55398871 1.5539886... 
3.51185 59486 3.51185 59497 3.51185 594 3.51185 594. 
6.961 53 30856 6.961 53 30856 6.961 53 308 6.961 53 308 
12.31107 74009 12.31107 74009 12.31107 73 12.31107 73 
20.19898 86459 20.19878 86459 20.19898 86 20.19898 86 
" 32.2288915016 32.22889 15016 32.2288915 32.2288915 


Now consider the calculation of the eigenvalues from computed values of 
det (A — AI) working with a mantissa of 30 bits. Since 2-* is approximately 
10-® we might expect that the larger eigenvalues would be obtained accurately 
using HyMan’s methcd of evaluation (or any other ‘‘good” method of determinant 
evaluation) but that the small eigenvalues would be hopelessly inaccurate. 
Moreover if we found this to be true we would not regard it as a criticism of the 
method. Such shortcomings are, in general inveitable. However we find that 
using HyMAN’s method of evaluation, all eigenvalues are obtained very accurately. 
The computed values are given in Table 6. 

On the other hand consider the matrix B, below. It obviously has the same 
eigenvalues as A,, since 

det (A,, — AI) = det (B,—Al). 


n n—1 
(n—1) (m—1) (m—2) 
p, =|" —2) (#—2) (n—2) (n—3) 
2 2 2 2 -'s 
ie.» 1 1 . t 








Yet if we apply HyMan’s method to the matrix B,, the small eigenvalues are 
determined completely inaccurately though, as may be seen from Table 6, the 
errors are in general only of the same order as those caused by changing the 
(12,1) element by 10°*. The result which requires an explanation is the first 
and not the second. 

The phenomenon is due to the special form of the matrix A,,. Although 
random changes of order 10-® in A,, alter completely some of the eigenvalues 
this will not be true of suitably correlated changes. This remark is of course 
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true of any ill-conditioned matrix but because of the special form of A,, the 
chance of correlated rounding errors is exceptionally high. 

Suppose we change all elements in the first column of A,, by the same quantity 
6,. It is obvious that the determinant is then (4+ 6,). If we change all the 
non-zero elements of the second column by the same 6, the determinant is 
unaltered and the same result holds for all later columns. Now in the evaluation 
of det (A — AJ) for any value of A, consider, for example, the computation which 
takes place on the 11th and 12th rows. For the eleventh row we calculate for 
some 4%, %2, %3,---, X44, flL[A Xx, +2 %.+3 XK X%y+---+(114—A) X x,,] and divide 
it by 11. For the 12th row we calculate 


PUTA x x2 X xo +3 X Hye +1444, + (12 — A) XH 9]. 


Up to the 10th term, the computations and the rounding errors involved are 
identical. We may therefore represent the sum up to the 10th term as the exact 
sum of 

1’ XX +2’ Xx%_+-+++10' X X49 


where 1’, 2’,...,10’ represent the same small perturbations of the elements 
1,2,...,10 in each case. Unless the contribution to the scalar product form 
the x,, term is very small compared with that from the first 10 terms, we can 
accommodate the rounding error made when adding the 11th terms by modi- 
fications to the elements 11, in each row. We can therefore predict that the 
calculated x; could have been obtained exactly, working with a perturbed version 
Aj of Aj, in which the perturbations are of the form 


jz 
3b, op o 

Oy oP os On 
6, 02 03 O4 
6, 02 03 4 


Since, as we have shown, perturbations in the positions (7,7+-1) are relatively 
unimportant, we see immediately that the perturbations are correlated in just 
such a way as to leave the value of the determinant substantially unaltered. 
Working with the matiix B,, no such correlation of rounding errors occurs. 

In [3] the evaluation of determinants of matrices obtained by rearranging 
(Ay.—AI) was discussed, including the evaluation of det(C,,—AJ) below, by 


P12—A M1 10 rar 
14 41—A 10 rey 7 
100—A 9 1 


(Cy, — 41) oa 


eo e2 8 oe &@ 0 8 6 @ 6 & DD 








Gaussian elimination with interchanges. There it was suggested that the loss 
of accuracy for small values of A may be due to the fact that interchanges take 
place at each stage after the first and as (page 31) a result the final pivotal 
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element is very small. The following analysis shows that interchanges play no 
important part in the loss of accuracy. 

We display below the first step in the evaluation of the determinant for 
A=0 and for economy of presentation we have used floating decimal arithmetic 
with 6 figures in the mantissa. The first step is the subtraction of .916667 times 
the first row from the second. Because of the Hessenberg form of the matrix, 
rows 1 and 3 to 12 are unaltered. The comnutation has been performed observing 
true floating-point rules so that, for example, the second element in the reduced 
row is given by 

f1(44 — 11 x.916667) = 11 — 10.0833 = .916700. 


The reduced row is 


0 .916700 .833300 .750000 .666660 .583330 .500000 .416670 .333330 
.250000 .166670 .833330x107. 


We now calculate exactly the effect of the rounding errors made up to the end 
of this single stage. The reduced matrix has rows 1 and 3 to 12 exactly as they 
were originally and row 2 is that given above. The determinant of this reduced 
matrix will be unaltered if we multiply the first row exactly by .916667 and 
then add it exactly to the second row. This gives a matrix identical to C,, apart 
from the second row which is exactly 


11.000004 11.000037 9.999970 9.000003 7.999996 6.999999 6.000002 
5.000005 3.999998 3.000001 2.000004 1.000000. 
Because of the simple structure of the matrix C,, the determinant of the perturbed 
matrix may be readily computed and we obtain 


det (Cj2) = 1+ 134 


the second of the terms being the contribution from the rounding errors. For 
computation with 8 and 9 decimals the rounding errors are exactly 45th and 
roo th those given above. Even with 9-figures in the mantissa the rounding 
errors in the first stage change the determinant completely. In fact if we con- 
tinue this analysis we find that the rounding errors made in the remaining 
10 stages of reduction have a smaller total effect than those in the first. We 
see then that the large error in the determinant should not be ascribed to the 
interchanges or to the cumulative effect of the rounding errors made in all 
11 stages. In fact no interchange is performed in the first step and this is easily 
the most damaging. Further it is precisely because the numbers in the successive 
derivatives of row 2 become progressively smaller that the rounding errors in 
the later steps are less important than those made in the first. 

On the other hand if we perform Gaussian elimination on the matrices 
(Dy. — AI), where D,, is defined by 


kt 1 7 
$9.4... 2 2 
29 . 
— 3 3 ae | 


.e «.53.¢0¢ £ € 2£ 4 2s Oe Ss OS 2 4 oe SS 
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the rounding errors are again correlated. For example if we take A=.123456 
and use floating-point arithmetic with 6 figures, an interchange takes place in 
the first step. The multiplier is .876544 and the reduced first row is 


0 —.644870 —.753090 —.753090 —.753090 ... —.753090. 


Proceeding as before, we see that the first reduced matrix is the exact equivalent 
of (D,, — .123456/) except for the first row which is 


876544 .99999987776 .999998 .999998 ... .999998 


and the diagonal elements, which have errors because we cannot subtract .123456 
from 2,3,...,12 exactly in the domain of six figure arithmetic. The small 
diagonal errors are of no importance since the determinant of D,, is not very 
sensitive to changes in these elements. The errors in the last ‘ten elements of 
first row are all exactly .000002 and are therefore correlated in a way which 
makes them unimportant. 

The above analysis shows why it is that the smaller eigenvalues, tound as 
zeros of det (Cj,— AJ) evaluated by Gaussian elimination with interchanges, are 
completely inaccurate while the same eigenvalues found as zeros of det (D,,— AJ), 
evaluated in the same way, are very accurate. 

Finally by simple row operations, as was shown in [3], we may reduce 
det (C,, — AI) to det Ky. where Ky, is the triple diagonal matrix shown below. 





r1i—A A 0 we ® 0 
44 1-A A a 0 
wid... 
ng Peep eer g oh oe tear 
21-A A 
a 4 1-A_ 





As in §3 we may show that the evaluation of the determinant of this triple 
diagonal matrix for any permissible value of 4 will give the exact determinant 
of a perturbed matrix which is also triple diagonal. We denote the perturbations 
by «,; on the diagonal, 8; on the super-diagonal and y,; on the subdiagonal; we 
may readily obtain very satisfactory bounds for «,, 8B; and y;. Now by reversing 
the row operations which were necessary to produce the triple-diagonal matrix 
from C we may find the perturbations of C,, which are equivalent to those 
in K,.. The matrix of these perturbations is easily seen to be Pj, given below. 


(04 + Ya) «++ (Mo +Biot+Yn) (%+Biut+ vie) (2 + Biz) | 
Ye ves (G10 + Bio t+ Yu) (411 +Biu+ Vie) (12 + Bra) 
dap eice.* + ees eh doth eat ce ti¢ 84 =a 
0 0... (0+ Y1) (41+ But i2) (%2+ Are) 
00 ... Yu (11 + 12) (12 + Biz) 
o > Owes 0 Viz M2 | 








Again we have a correlation of errors which nullifies their effect and eigenvalues 
found from the triple diagonal form are very accurate. 
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Although examples of this type provide interesting exercises for error analysis 
the results are irrelevant to the assessment of a numerical technique. The methods 
which gave accurate results for the smaller eigenvalues are no ‘‘better’’, in general, 
than those which did not. 

Acknowledgements. The work described here has been carried out as part of the 
Research Programme of the National Physical Laboratory and this paper is published 


by permission of the Director of the Laboratory. The author wishes to thank Dr. 
E. T. Goopwin for reading the manuscript and making numerous valuable suggestions. 


(1] 
[2] 
[3] 


[4] 
[5] 


[6] 


[7] 
[8] 


[9] 
[10] 
(11] 


(12) 


References 


Carr, J. W.III.: Error Analysis in Floating Point Arithmetic. Commun. of 
the Assoc. for Comp. Mach. 2, 10 (May 1959). 

FiscHER, P.C.: Automatic propagated and round-off error analysis. (To be 
published.) 

FRANK, W. L.: Computing eigenvalues of complex matrices by determinant 
evaluation and by the methods of DANILEWskKI and WIELANDT. J. Soc. 
Indust. Appl. Math. 6, 378 (1958). 

Givens, W.: Numerical computation of the characteristic values of a real 
symmetric matrix. Oak-Ridge National Laboratory, ORNL-1574. 

Hyman, M.: Eigenvalues and eigenvectors of general matrices. Presented at 
the 12th National Meeting of the Association for Computing Machinery, 
June 1957, Houston, Texas. 

Lipski, V. B.: O sobstvenyh znacéeniyah summy i proizveniya simmetriceskih 
matric. Doklady Akad. Nauk 75. 769 (1950). 

METROPOLIS, N.: Automatic round-off error analysis. (Unpublished.) 

White, P. A.: The computation of eigenvalues and eigenvectors of a matrix. 
J. Soc. Industr. Appl. Math. 6, 393 (1958). 

WILKINSON, J. H.: The calculation of the eigenvectors of codiagonal matrices. 
Computer J. 1, 90 (1958). 

WILKINSON, J. H.: Rounding errors in algebraic processes. Proceedings of Inter- 
national Conference on Information Processing. Unesco 1959. 

WILKINSON, J. H.: The evaluation of the zeros of ill-conditioned polynomials. 
Numer. Math. 1, 150. 

WILKINSON, J. H.: HousEHOLDER’s method for the solution of the algebraic 
eigenproblem. Computer J. 3, 23 (1960). 


Nationai Physical Laboratory 
Teddington, Middlesex/Great Britain 


(Received June 1, 1960) 








— 


Numerische Mathematik 2, 341 —343 (1960) 


Bestimmung der Eigenwerte 
des allgemeinen Eigenwertproblems 
Von 


HANS-LOTHAR HAAS 


1. Aufgabenstellung 
Das allgemeine Eigenwertproblem 


{AB+C}r=4 (1) 
14Bt sich dann auf das spezielle Eigenwertproblem 
E+ Wr=q (2) 


leicht transformieren, wenn % regular ist. Es wird lediglich Gleichung (1) von 
links mit 8 multipliziert. Fiir singulare Matrizen 8 soll im folgenden ein Vei- 
fahren gezeigt werden, das die Transformation auf das spezielle Eigenwert- 
problem bewirkt. 
2. Transformationen 

Nach EGERvARY [J] kann jede Matrix als Produkt zweier Matrizen dargestellt 
werden, die beide den Rang der urspriinglichen Matrix besitzen, deren einer 
Faktor eine Dreiecksmatrix, und deren anderer Faktor eine Dreiecksmatrix oder 
ein Teil einer Dreiecksmatrix ist, je nachdem ob die urspriingliche Matrix qua- 
dratisch oder rechteckig ist. In einer modifizierten Art fiihrt die Zerlegung auf 
das Produkt einer regularen unteren Dreiecksmatrix, der zugehérigen Rang- 
einheitsmatrix und einer regularen oberen Dreicksmatrix. 


V=U-E,-D (3) 

U: Untere Dreiecksmatrix 

®D: Obere Dreiecksmatrix 

€,: Rangeinheitsmatrix. 
Eine der beiden Dreiecksmatrizen ist im allgemeinen ungeordnet, desgleichen 
die Rangeinheitsmatrix. Die reguliren Dreiecksmatrizen, deren Rang eventuell 
gréBer ist als der von % und von denen in diesem Fall gewisse Zeilen oder Spalten 
nicht eindeutig bestimmt werden kénnen, ergeben sich aus den singularen Dreiecks- 
matrizen am einfachsten durch Ersetzen der verschwindenden Hauptdiagonal- 


elemente der geordneten singularen Dreiecksmatrizen durch die entsprechenden 
Elemente der Einheitsmatrix. Die geordnete Rangeinheitsmatrix hat die Gestalt 


€ O 
€, = 

hr 
(©: Einheitsmatrix 
©: Nullmatrix, 
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wobei eventuell die letzte Zeile oder Spalte verschwinden kann. Die Rang- 
einheitsmatrix hat stets den gleichen Rang und die gleiche Zeilen- und Spalten- 
zahl wie die zugehérige Matrix %. 

2.1 Transformation von 8 auf die Rangeinheitsmatrix ©, , 

Gleichung (1) wird so transformiert, daB 8 in die Rangeinheitsmatrix &, , 
iibergeht. UAAB+CGO}D +z, =q 

{AG +C}n =q. (4) 

Es werden nun zwei Fille unterschieden 

1. G, , hat den Héchstrang (7, =n) 

2. &,, hat nicht den Héchstrang (0<7,<) 

r,: Rang von $ und &,, 
n: Zeilenzahl = Spaltenzahl von 8, G, &,,. 

Der dritte Fall 7,=0 tritt nicht auf, weil dann das Problem von vornherein 
als entartet ausgeschieden werden kann. 

2.1.4 Lésung fiir 7,=n 

€,, hat den Héchstrang und damit ist die Transformation auf das spezielle 
Eigenwertproblem bereits beendet. Dic Eigenwerte und Eigenvektoren kénnen 
nach bekannten Methoden bestimmt werden. 

2.2 Transformation fiir 7,;<" 

Die Koeffizienten von Gleichung (4) werden in Blockmatrizen zerlegt, wobei 
angenommen wird, daB die Rangeinheitsmatrix €, , geordnet sei. 











Gi) G q 
A =]|'|. 
tale) + lelta=(" (5) 
Gleichung (5) wird so transformiert, daB ©; in die Rangeinheitsmatrix €,, 
iibergeht. (cs 7 1) q] 
waaay”? 1} D2, = 
| ja} * cl aa 
Bi’ cy’ a] 
(S-Ele-f 
| - E,» : q 











Auch hier sind zwei Falle zu unterscheiden 
1.6, hat den Hochstrang (7,=" — 7) 
2. €,, hat nicht den Héchstrang (7,<" — 1”) 
r,: Rang von G3 und &, 4. 
2.2.1 Transformation fiir 7,=n —1, 


Aus dem unteren Teil der Gleichung (6) folgt, daB ein Teil der Variablen 
verschwindet. Mittels einer weiteren Transformation bzw. durch Streichen der 
letzten 7, Zeilen und der zu den verschwindenden Variablen gehérigen Spalten 
geht Gleichung (6) iiber in 


{A Na +" 313=4q, (7) 


die wieder die Gestalt von Gleichung (1) hat, deren Koeffizientenmatrizen aber 
statt m Zeilen und Spalten nur mehr —y, Zeilen und Spalten besitzen. Um 
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das spezielle Eigenwertproblem zu erhalten, wird Gleichung (7) weiter nach 
Abschnitt 2.1 ff. transformiert, bis entweder €,, den Héchstrang oder €,, nicht 
den Hochstrang hat. 

2.2.2 Lésung fiir 7,.<"—1”. 

Die zu ry gehdrige Koeffizientenmatrix der Gleichung (6) enthalt mindestens 
eine verschwindende Zeile und daraus folgt, daB die Lésung mindestens einen 
freien Parameter hat, d.h. jedes A ist Eigenwert und jedes rx, Eigenvektor. In 
diesem Fall hat das Problem keine eindeutige Lésung. 


3. Eigenvektoren 
Die Eigenvektoren des allgemeinen Eigenwertproblems ergeben sich aus den 
Eigenvektoren des zugehérigen speziellen Eigenwertproblems durch Riicktrans- 
formation. Infolge der dimensionsvermindernden Transformationen bei der 
Reduktion des allgemeinen auf das spezielle Eigenwertproblem hat das allgemeine 
Eigenwertproblem keine linear unabhangigen Eigenvektoren mehr, sondern 
nur mehr soviele wie das zugeordnete spezielle Eigenwertproblem. 


4. Verringerung des Rechenaufwandes 

Eine Verringerung des Rechenaufwandes ergibt sich, wenn die Transforma- 
tionen nicht bis zur Rangeinheitsmatrix sondern nur bis zur Dreiecksmatrix 
durchgefiihrt werden, da schon an dieser Stelle die entscheidende Rangfrage 
geklart ist. Des weiteren kénnen die Matrizenmultiplikationen und -inversionen 
umgangen werden, wenn statt der Transformationen der GauBsche Eliminations- 
algorithmus gleichzeitig auf beide Matrizen angewendet wird, aber nur eine der 
beiden Matrizen auf Dreiecksform gebracht wird. 


5. SchluBbemerkungen 

Fiir die Berechnung der Eigenwerte mittels eines Digitalrechners eignet sich 
die Gau8-Elimination am besten, weil dabei die Anzahl der notwendigen Opera- 
tionen am geringsten ist. AuBerdem treten unabhangig vom gewahlten Verfahren 
bei der Rangbestimmung infolge der unvermeidlichen Rundungsfehler die be- 
kannten Schwierigkeiten auf. 

Zusammenfassung 

Es wird ein Verfahren hergeleitet, das mittels rangerhaltenden aber eventuell 
dimensionsvermindernden Transformationen das allgemeine Eigenwertproblem 
auf das spezielle Eigenwertproblem transformiert. Das allgemeine Eigenwert- 
problem darf dabei singulare Matrizen enthalten. Hierbei kann der Fall ein- 
treten, daB das Problem keine eindeutige Lésung hat. 
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On the Kantorovich Inequality 
By 
ANDREAS H. SCHOPF t* 


Let A be a positive definite hermitian matrix whose characteristic roots are 


(1) Ay 2 d4g2-:-24,>0, 

let 

(2) x2 = A,/A, 

be the spectral condition number, and let x=-0 be any nonnull vector. Define 
(3) HM, = Ml, (x) = x¥ A’ x, 

where the superscript H represents conjugate transpose. Then 

(4) 1S ty gr fy alp S (xe +974. 


The first inequality is well known and can be proved in many ways. The second 
is one formulation of the Kantorovich inequality [2]. Several proofs are to be 
found in the literature. The one given by Kantorovich himself makes use of a 
spectral decomposition and the solution of a minimizing problem. More recently, 
GREUB and RHEINBOLDT [J] have given a proof applicable in a more general 
normed space. 

Since in the theorem there enters explicitly only a particular two-space, 
whatever the dimensionality of the entire space, it is reasonable to expect that 
a proof might exist that never requires explicit consideration of more dimensions. 
The proof to be given here is of that character. Before giving it, however, 
certain consequences of the theorem itself may be of interest. 


Let X be a matrix of any number of linearly independent columns, and let 
(5) M, = X# A’X. 
By a natural generalization, define 
(6) u,(X) = det M,. 


In fact, if X is of rank o, and if B®) designates the o-th compound of the matrix B, 
then 
(7) u,(X) = Me) = X@H Alor Xt) 


and X) is a vector. Hence yu,(X) is the same as the moment yu, (X)) computed 
with the matrix A‘) in place of A. Since the characteristic roots of A® are the 


* Deceased, September 29, 1959. The actual writing and organization of this 
paper was done by the undersigned, making use of scattered notes, and of such 
fragments of conversation as he can recall. A. S. HOUSEHOLDER. 





On the Kantorovich Inequality 345 


products o at a time of those of A, the condition number is 


(8) 9 = Ay.» Agl(An—ota>++ Ag) 
and the inequality becomes 
(9) 1S pty 4 (X) fy a (X)|p2 (X) S (% + 3)7/4. 


In particular, let y=0, and take the columns of X to be any @ columns of J. 
Then 1, (X) is a principal minor determinant of A of order @, u_, is the correspond- 
ing principal minor determinant of A, and uj»=1. Hence, if the theorem of 
corresponding minors is applied, the result is 


(10) 4 A(1,...,0) A(o+1,..., m)/A(1,..., 2) S (x, +%9)7/4, 
where A(i,7,...) is the determinant of the principal minor formed from the 
i-th, j-th, ..., rows and columns. The inequality on the left is again well known; 


that on the right was suggested by Ky FAN (private communication). 
Primary interest lies in the inequality on the right of (4), although one way 
to establish that on the left is to observe that the matrix 


G A’-1 (A x, x) es (* by ) 
My My-1 


% 
is at least semidefinite, and is definite unless x is a proper vector. Hence if 


(41) Vy (%) = M,(%)/My—1(*), 
then 

Yv+1(%) = V(X) 
Furthermore, 
(12) A, 2 ¥,(x%) 2A). 


Now consider in particular 
Ve = (A) =7,(4¥ +444). 


Since 
(Mya + 2A pm, + A? oy 41) Yy — (My + 2A M41 + AP My 42) = 9, 


extreme values of y, are the roots of the quadratic 
det ae My41— VM = 0. 
Poti — VP» Py+2— Vy +1 
Let these roots be y’ and y”’. The function 
p(a,B)=4aB/(a+B)?, aB>O, 


(13) 


satisfies 
YP (a, a) =1, 


and decreases monotonically as «/B>1 increases. But 


4y¥(%)_ _[Y¥v4+1(*) —Yo(*)) [¥r-42(*%) — Yr t1 (¥)] 





PO) = ale) ye+a(*)—r(2)]2 
= 9[Yr42(%) — Yr41(%), Yr4a(%) — Yy(*)] Yy(%)/Yy42(*) 
= Vy» (%)/Yo4a (x) , 
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where use has been made of the formula for the sum and product of the roots 
of the quadratic (13), followed by (7). However, in view of the monotonicity 
property, (12) implies that 


p(y’, ¥") = P(A, A,)- 


Y» (x) /Yr41 (x) = P(A, A,) ’ 


and reciprocation yields the right member of (4). This completes the proof. 


Hence 
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A numerical study of Webet’s real 


class number calculation 
Part I 


By 
HARVEY COHN 


1. Historical background and summary of results 


The calculation which concerns us here arose in HEINRICH WEBER’S celebrated 
proof of KRONECKER’s theorem that each abelian field over the rationals is a 
sub-field of a field generated by some root of unity [7; pp. 762—783]. 

In the course of this proof, one major obstacle which WEBER overcame was 
a lemma that the field, of degree 2‘~', 


(1.1) R, = R(exp 22/2'), t>2 

has odd class number. This, in turn, led to the study of the real sub-field of 
degree 2'~? 

(1.2) f, = R(2cos 22/2'), t>}3. 

WEBER developed class number formulas which have been scarcely improved. 
The class number of f, is given by 

(1.3) h, = 83 84--- 81, 


where g, is defined as the imdex of the special multiplicative sub-group of real 
units of f, generated by the values 


(1.4) T, = tan 5° 2/2’ 


within the set of all normal units of f, (as defined below, § 2). We concern our- 
selves entirely with f, and h,. 

We might note, however, that the class-number of &, is given by a formula 
of the type h, k, where the calculation of k, is relatively straight-forward. The 
current results are [2] 


kg = 1, (Classical) 
(1.5) 4 ky=hg=1, he=17, hp =17X21,121 (WEBER) 
kg = 17 X 21,121 X29, 102, 880, 226, 241 (Hasse) 


This would be an appropriate point to remark that these values of k,, how- 
ever large, are not as suitable a problem for large scale computers as the values 
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of h,. In the case of the k, the calculation is rather straightforward (although 
it involves a formula [2] too long to reproduce here). The calculation of the k, 
would involve multiple precision difficulties more than any other kind for the 
programmer, and the answers known even now are too large to be appreciated. 


The calculation of h,, in contrast, involves the index of the quotient of two 
infinite groups (of units, see § 2). The problem can be “‘visualized”’ geometrically 
as the index of a sub-lattice in a lattice in 2'~* dimensional space, but the question 
of numerical accuracy would seem formidable since conjugates of a unit necessarily 
tend to remain ‘‘unscalable’’. The type of numerical analysis we perform is the 
well-known method [J, 5] of exact integral arithmetic modulo # for special values 
of the prime p. This is in the theoretical sense an ‘‘approximation modulo p”’ 
in the same way that decimal arithmetic is an “‘approximation modulo oo”. 
Thus the calculation of h, requires a large scale computer all the more for being 
conducted in integral arithmetic since # ranges up to a few million and there is 
an extensive use of indices and primitive roots. 

Very little is known about the f,. For the field ®,=R(2#), easily hs=1. 
WEBER proved /,=1, from REUSCHLE’s tables of ideal factors rather than by 
his own method (1.3). We briefly consider WEBER’s argument [7; p. 808]: 
Letting the discriminant of f, be D,, then by MinKowskI’s theorem each ideal 
class contains an ideal of norm < N, where [7; pp. 690, 756], 


a om 
| N= er 


(1.6) 


Thus N,=4.24..., N;=111.36..., Ng=881,830. ..., etc. Now it is an elementary 
verification that in f, the ideals of norm <4 are (1), (2 cos 27/2*), or (2 cos 27/24)? 
which are principal ideals of norm 1, 2, or 4 (norm 3 being excluded). For f,, 
indeed, REUSCHLE’s tables [4; pp. 451—2] show all ideals of norm <=1000 to 
be principal. Thus REUSCHLE’s tables have sufficient data to establish that 
h,=1. (This fact seems to have been previously overlooked. Indeed, WEBER 
[7; p. 808] thought that REUSCHLE had disposed of only the ideals of norms 
<= 100.) 

The case of f, is clearly hopeless to approach, similarly, by finding ‘‘enough”’ 
unique factorizations to satisfy MINKOWSKI’s theorem. 

Now WEBER casually conjectured [7; p. 808] on the basis of kg=17 that 
hg>1. The evidence in REUSCHLE’s tables does not substantiate this inference 
at all [4; pp. 455—458] and, indeed, in the present calculation we incidentally 
establish that h, is not divisible by 17. 

Actually, we shall show: 


he=1, or 16015 h, = prime S 83,921; 
h,=1, or 1601Sh,; 

h,=1, or 1409S h,; 

h,=1, or 257Sh,; for t=>9. 


We know h,=1 for #5 from before. 
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All these results fail to contradict ‘‘h,=1’’. Yet we might believe that ulti- 
mately h,>1 (since we have a chain of fields f;¢fy¢ ---). Thus we might be 
prepared to believe that the much-sought case, where the field f,= R (2 cos 27/2') 
finally does have a class number >1, might involve a fantastic class number, 
perhaps larger than the ‘‘complex”’ contributions k,, that were just enumerated. 


2. Summary of relevant portions of algebraic number theory 


We summarize well-known results by first introducing the convenient symbols: 


(2.4) m=2' (> 32), 
(2.2) n= m/4 = 2-8, 

(2.3) l= m/8 = 2-*. 

Let 

(2.4) ¢ =exp2nai/m=Cp, 
(2.5) s= 00. 

Then 

(2.6) Csrr=—s, 

(2.7) c=, 


~ 


which is a simple consequence of the ‘‘rational’”’ results that 


(2.8) 5"=1 (modm), 
(2.9) 5"=1 (mod2m). 


We define the Kummer normal units as in (1.4), by 





(2.10) t= 41> 2 tt = tan §*x/m, m0, 8)...,8—4: 
2S 

then 

(2.11) T5414 = 1/T,. 


Now the field R(¢+ 1/C) = R(2cos 22/m) is of degree m and an arbitrary algebraic 
integer is an arbitrary polynomial in (¢+1/¢) with integral coefficients. The » 
conjugates are generated by the substitutions (2.5) ¢,+¢,, where s=O,1,..., 
(7 —1). These substitutions pair off according to (2.6), since /=n/2. 


We define a normal unit as a unit / (2) in f, for which, (as in (2.11)), 


WEBER’s famous lemma [7; p. 818] makes the factor g, (of (1.3)) equal to the 
index of KUMMER’s normal units as a sub-group of all normal units. 
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Hence, by general results in linear dependence, the prime g divides g, if 
and only if integers a), a,,..., @;_, exists, not all divisible by g, such that 


(2.13) TT +++ Tt = YG 


for some algebraic integer y=yy=y (Co), (in the usual notation for conjugates 
in f,). 

The theory becomes differently embellished for the purpose of other calcu- 
lations. (We cite references [1] and [5] for a variety of results.) 


3. Use of modular arithmetic 
We next consider an odd prime g which is to be a hypothetical divisor of g,. 
For this prime g we consider another prime # such that 
(3.1) p=1 modmgq. 
We let g be a primitive root of p. Now pis a prime of first degree in f,, [7; p. 758]; 
thus 
(3.2) P= PoPi--+ Pai: 


We then know that every algebraic integer is congruent to m (distinct) rationa 
integers modulo p,, as s=0,1,...,("—1). In particular for some conjugate 
say s=0O, 


333) g* =2=C (mody,). 

We next define the rational integers T,, by 

(3.4) it,=T,=—1+2/1+2" (mod pp) 
and the indices L,, modulo(p — 1), by 

(3-5) g’*=T, (mod). 


We now proceed to write down the conjugates of equation (2.13). To make 
the notation simpler, with no real loss of generality, we shall take m= 32, 1/=4, 
then since T,,4=— 1/t, we have 


wt wr a =r 
We ty wt =— yf 
wy wit = 


—a —a. —fe .... 
3° TT, * te * = — yf. 


3.6) 


Now, if we take the system modulo p, and note 
(3-7) (—)'y,=G,=g" (mod p,), 


we find that the following type of congruence emerges from (3.6) modulo py, 
and, hence mod # since all terms are rational: 


(3.8) {Ty Te: Ty Ty} = G4" (mod), ete. 
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Here we took the fourth power to remove the factor 7 arising through (3.4). 
Then taking indices in (3.8) and like equations, we find 


Qo l,+4,L,+4,L,+ 4,L3= 9X, mod(p — 1)/4 
ao Ll, + a,L,+a,L,;—a,L,= 9X, mod(p — 1)/4 
a) L,+ a, L3—a,L,—a,L,=9qX, mod(p — 1)/4 
ay Ls — a, L,—a,L,—a,L,= 9X, mod(p — 1)/4. 


(3-9) 


We note that ®, the determinant on the left, is ‘‘almost’’ a circulant. By the 
usual technique, we can show it to have the factors (Lp+L,&+L, &#+L, &) 
where & takes on the four primitive eighth roots of unity. This method generalizes 


easily to any m. 
More generally, g can be a divisor of g, only if 


(3.10) q|® 

where, for € a primitive n™ root of unity, 
j-1 J-1 

(3.44) +£O= fT UL,ee. 
h=0 s=0 


We note, however, that + @ is the norm of 
(3.12) Po=Llot+hét+-:-+L,1,8" 


over the field of the »™ root of unity. Hence by general theorem of ideal fac- 
torization, g|® if and only if either (or both) of the following are valid: 


(3.13) L=l,=::-=L1,1=0 (modq), 
(3.14) g=1 modn. 


We exclude the first alternative (see § 4), and note that equation (3.14), like 
equation (3.1), leads to: 


(3.15) Y= %1-+-U-1 

(3.16) Es x= i" (mod aa), 
where 7 is a primitive root of g. Thus if we set 

(3.17) f(x) =LotLyx+--+L,,x, 
(3.18) h=t(""), 

(3.19) F=foh,---h-1 


we find that if the prime q satisfies equation (3.14) then q can be a divisor of g, 
only tf 
(3.20) q|F. 


The problem is now wholly in rational arithmetic. 
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4. The excluded cases 
The theoretical scope of the present paper (Part I) is enhanced by the follow- 
ing result: 
Lemma 4.1. Any odd prime q which divides g, must either satisfy condition 
(3.14), g=1 mod 4, or else 


(4.4) bed 


for KUMMER’S normal unit tT. 
To see this, again take m= 32, for convenience, and note the equations (3.6) 
with the t, kept in order (rather than the a,): 


wont my pa 
TY mm m= — 
(4.2) i Sie Si Ba 
T *T, “Ty T= 3 
To “1 Ty * T, “* T° = — 7. 


Now, A, the determinant of the a, is very much like the determinant ® of the 
L, in system (3.9). Its factors are conjugates of (a3+ a, &+ a, &*+ ay &). Next, 
in (4.2), we concentrate on (say) the first column and raise each line to a power 
equal to the corresponding cofactor A;. Then on multiplying all lines together, 
we find 


(4.3) T= +(e vt vets)! 


Finally, as in § 3, if g =1(mod m), then gtA, since not all a, are =0 (mod q), 
leaving as the only alternative (4.1). Q.E.D. 

These remarks fit in with our use of the computer for purposes of excluding 
possible divisors g of g,. Hence if we restrict ourselves to values of g satisfying 
(3.14), then the only cases we miss are those where T, is a perfect g power. The 
excluded cases are more easily subject to decimal arithmetic, which we shall 
pursue in Part II of the present paper. We might note, for the time being, 
that the possibility of (4.1) might even be always excluded by a theoretical 
argument (which unfortunately seems too remote at present). In any case if 
t=y! the index desired becomes divisible by g’ which is =1modu by the 
Fermat-Euler Theorem. Thus, whatever happens, 


(4.4) g,=1 modu. 


In one way, very little harm is done by the excluded cases if we consider 
the sheer size of the factors g,. For if m=64, then /=8, and hence the minimum 
contribution from a perfect power is q’=3*=6,561, far beyond the range of q 
that we have been able to reach in this part of the paper. The restriction (3.14) 
will, of course, enormously cut the eligible g and speed the calculation. 


5. The machine calculation 
We first define the function 


(5.1) M(a) = 2*||(a —1), 





4 
| 
! 
* 
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where 2*| (a — 1) while 2*** +(a—1). Then we consider sequentially all odd primes 
(5.2) g=1 (mod8). 


These primes, considering our exclusion, are eligible to divide g, if and only if 
2'|4M(q). This leads to a set of eligible 


(5.3) m = 32, 64,...,4M(q). 
For each m and q we find the first odd prime. 

(5.4) p=1 (modmg). 

Of course the same # could suffice for several m, since 
(5.5) M(p) =m 


from (5.4), and often M(p)=2m, or 4m. In practice, for instance when g= 769, 
p= 590,593 for m= 64, 128, 256; (see Table 1). 

Thus the imput data consist of the primes g listed with the smallest ; satisfy- 
ing (5.4) for m; from 32 to m,=4M(q): 


q, Pi, (= 32), 
q, Pe, Me(= 64), 
(5.6) ) 7) Ps, M3 (= 128), 


e+ * Oe Ow 6 





The imput data fill columns 1,4, and 7 of Table 1; column 2 determines the 
number of m, that are required. The imput data were calculated by machine 
to ensure accuracy. 

The progress of the machine calculation was as described in § 3: The machine 
calculates the primitive root g (see §6 below); then z is calculated in (3.3) 
modulo #, actually, and with it the set of / values O<s<l, in reduced positive 
residue classes: 


(5.7) z,= 2" = (z,_,)° (mod). 

As a check, reducing all residues to least positive residues, 
(5.8) 4+%=p. 

Then we calculate and store 

(5.9) K,=2z,—1. 

(5.40) M,=2z,+1. 


Next the machine generate the powers of g and, by a ‘“‘table look-up” finds 
the index of each K,, and M, (mod #—1). Then we store L, found by sub- 
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tracting these indices (mod # — 1) ; 

(5.44) 0OSL,Sp—41, 

(5.42) g’* = K,/M,=T, (mod). 


Finally, turning its attention to g, the machine calculates the primitive root 7 
and x as in (3.16) as well as the /, of (3.18), all modulo q. 


The main “‘yes or no” problem is whether q|F, as in (3.20). This test is also 
internal. In either case, the machine produces as output (in decimal): 


pPeqi 
m z x (blank) 


(The / values of K;,) 
(5.43) (The / values of M,) 
(The / values of L,) 
(The / values of /,) 
(Signal to indicate p| F or gtF). 


The problem is, incidentally, quite amenable to the four-column output. 


6. The primitive root sub-routine 


The primitive roots are found by trial and error. The machine generates 
p—1 (or g—1) powers of a trial value g and verifies that, (restricting the de- 
scription to ), 


(6.1) g*=1 (mod) 


exactly when x=~—1. (This also serves as a check for machine error.) Of 
course, not every g is used. By the theory of quadratic residues we ignore g=2, 
when #=1(mod 8). We also ignore any other g for which 4g|(p—1) as well as 
any g divisible by a square. 

Thus, looking over our results, when $= 56,929, where g=14 (a rather large 
value), we had only to first try g=5, 7,10, 11,13. Actually relatively little 
time is spent on these trials, so that in retrospect, there was no need to encode 
the usual techniques for building a primitive root as the product of non-primitive 
roots. Indeed, we increment g from 2, yielding the Jeast primitive root. 


The sub-routine has a switch for making a record of z (or x) while the powers 
of g (or j) are enumerated, according to (3.3) (or (3.16)). 


VINOGRADOV’s table of least primitive roots [6] was used for spot-checks. 
Interestingly, one error in the table emerged: for /=3631, the first primitive 
root g is 15. The value g=21, listed in this table, was seen to be the second 
primitive root. (Undoubtedly this table was computed by hand!) 

















A numerical study of WEBER’s real class number calculation 355 


7. Checks and tests 


Certain checks were built into the program (after the primitive root sub- 
routine). Among these were (5.8), and a double system of checks for the look-up 
of indices of K; and M;. Here powers of g were generated all over again and 
a check was made on the number of “‘recognitions’’ or the number of times 
g*=K, or M,. The total must, of course, be 2/=mn. Even so, after the last 
recognition, powers of g were continually generated to ensure that equation 
(6.1) is valid precisely at x=) —1. 

The programming had very simple loops and was no great logical problem. 
The difficulty lay in the need for very tight loops in the look-up of K, and M, 
(which requires practically all the computer time). To make testing difficult, 
Jacosi’s tables of indices [3] do not go beyond #=1000, while our first case 
is p=5,441. As a result, two test calculations were made with imputs: 


q=5, p=2441, m=16, (n=4), 
q=29, p=920, m=32, (n=8). 


The first is legitimate, but the second violates (4.2). In either case the problem 
was permitted to run and the memory dump was checked against hand calcula- 
tions made on the basis of JAcoBI’s table. 


8. Remarks on the running of the problem 


The problem was run on GEORGE, with a memory of 4096 words of 40 bits 
of magnetic core storage. The word structure is two-address and sequence- 
controlled with amazing versatility. Among other things, an actual address can 
be used as effectively as a memory location of the address. The speed is roughly 
50,000 two-address operations per second. 


The running time of the cases was quite formidable, ultimately, since it 
approximated 5 #m/10® hours per case. Because of this fact, the problem was 
not feasible for (say) P>3-10® or m= 512. 


As a convenience in imput, where frequently two different m belong to one 
p and gq, we allow for a double case to be run in the memory. Thus two cases, 
which have the same # and gq, are subjected to the time-consuming look-up of 
L, at the same time. This economy involved only a minor change in the program. 
In two cases (somewhat whimsically), the smallest p was passed over in favor of 
the next # in order to use just one primitive root and one look-up for both of 
the cases. They are noted in Table 1, namely: 


q= 577, m= 128, p = 1,920,257 instead of 1,698,689; 
q = 673, m= 64, ~ =1,206,017 instead of 904,513. 


(See the * designation in Table 1.) 


The imput was on paper tape with program on one reader and data on the 
other. The output was also on paper tape. 
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Table 1. Main run 
1 2 3 4 5 6 7 8 9 10 
q | Mg | i p Mi) | ¢ | ™ z x Notes 
17 | 16 3 5441 | 64] 3 32 631 9 q|F 
5441 64 3 64 1638 3 
41 8 6 7873 64 5 32 516 27 
73 8 5 4673 64 3 32 2369 10 
89 8 3 11393 128 3 32 7007 37 
97 32 5 43457 64 3 32 13016 64 
43457 64 3 64 21195 8 
62081 128 6 128 22287 28 
113 16 3 3617 32 3 32 1084 18 
36161 64 3 64 19924 40 
137 8 5 30689 32 3 32 27 338 127 
193 | 64 5 30881 32 3 32 14966 43 
37057 64 5 64 28677 64 
49409 256 3 128 23 627 185 
49409 | 256 3 | 256 14096 125 
233 8 3 7457 32 3 32 2542 221 
241 16 7 38 561 32 3 32 19042 30 
46273 64 5 64 2150 111 
257 | 256 3 74017 32 7 32 72055 64 
82241 64 3 64 35841 249 
98689 | 128 7 | 128 81075 136 
328 961 256 3 256 230686 81 
1579009 | 2048 | ... 512 a a 0,06 
1579009 | 2048 | ... | 1024 
281 8 3 35969 128 3 32 33 568 60 
313 8 10 120193 128 5 32 6792 5 
337 | 16 10 21 569 64 3 32 15746 85 
21 569 64 3 64 16944 191 
353 | 32 3 33889 32 | 13 32 4458 283 
67777 64 5 64 40319 304 
858497 128 3 128 268 305 294 
401 16 3 51329 128 3 32 40059 45 
51329 128 3 64 23 288 268 
409 8 31 26177 64 3 32 18757 21 
433 | 16 5 83137 64 5 32 7 562 354 
83137 64 5 64 10785 238 
449 | 64 3 14369 32 3 32 5 566 122 
86209 64 11 64 38953 349 
517249 128 11 128 418330 221 
2643713 | 256| 3 | 256| 1828604 | 301 q|F 
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Table 1 (Continued) 
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1 2 3 4 5 6 7 8 9 10 
q M(q) i b M(p) | @& m | z | x | Notes 
457} 8 | 43 73121 | 32} 3 | 32 67463 | 207 | 
saa} 8 | 3 16673 | 32| 3 | 32 4368 | 345 | 
569 8 3 345953 32 3 32 278876 277 
577 | 64 5 36929 64 3 32 6 362 391 
36929 | 64| 31] 64 24173 | 171 
(1698689) | 128 3 128 | 1247461 400 next p 
*1920257 | 256/ 3 | 256| 492042 | 557 | used 
593 | 16 3 56929 | 32/| 14 | 32} 38344 | 392 
303 617 512 3 | 64 160953 384 
61 | 8 | 7 57697 | 32] 5 | 32/ 35828 | $509 | 
617| 8 3 59233 | 32) 5 | 32| 20407 | 478 | 
641/128 | 3 184609 | 32] 11 | 32] 112218 | 318 | 
697409 | 64; 3 | 64} 594875 | 601 
820481 | 256) 6 | 128| 720811 | 160 
820481 | 256| 6 | 256| 154505 77 
1312769 |2048 |... | 512 a mi 
673 | 32 5 64609 | 32) 35 | 32 22760 | 609 
(904513) | 64 3 | 64 56291 | 464 | next p 
*1206017 | 256 3 128 373844 118 | used 
761| 8 | 6 219169 | 32| 7 32| 180113 | 62 | 
769 | 256 | 11 270689 | 32) 3 | 32| 90931 | 40 | 
590593 | 256| 15 | 64) 147657 | 136 
590593 | 256) 15 | 128| 554154 | 304 
590 593 256 | 15 | 256 147670 | 85 | 
1968641 | 512| ... | 542 coe | oes | 
7874561 | 2048 | 1024 | 
g0o9| 8 | 3 25889 | 32| 3 | 32 5214 | 234 | 
857| 8 3 191969 | 32/ 3 | 32] 122002 | 669 | 
881 | 16 3 281921 | 64| 31] 32 69345 | 662 
281921 64 3 64 254561 | 767 
929 | 32 | 3 118913 | 128| 3 | 32 28146 | 671 
118913 | 128| 3 | 64 58144 | 889 | 
118913 | 128 | 3 | 128 51740 | 701 | 
937/ 8 | 5 149921 | 32| 6 | 32| 97210 | 67 | 
953 | 8 3 30497 | 32] 3 | 32| 17847 | 797 | 
977| 16 | 3 218849 | 32} 6 | 32| 96549 | 439 | 
375169 | 128| 13 | 64) 111219 | 629 | 
———EE7 —_ — EEE — - Zz —— | _ 
1009 | 16 11 64577 64 4-8 21402 762 
64577 64 3 | 64| 37246 179 | 
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Table 1 (Continued) 
1 2 3 4 5 6 7 8 9 10 
q | M(q) 7 p M(p) g m | z x Notes 
1033 8 5 198 337 64 5 32 75326 802 
1049 | 8 3 33569 |- 32 3 32 18449 | 461 
1097 | 8 3 140417 | 128 3 32 38745 486 
1129 8 11 325153 32 5 32 290042 692 
1153 | 128 5 405857 32 5 32 266162 1030 
442753 128 10 64 293 510 1086 
442753 | 128| 10 | 128| 233912 286 
2066177 | 256 3 | 256| 949550 943 
12397057 | 512 512 on Kas 
1193 | 8 3 38177 32 3 32 30679 524 
1201 16 11 192161 32 3 32 70958 343 
538049 | 64 3 64 255376 | 473 
1217 | 64 3 116833 32 5 32 85410 387 
155777 | 128 3 64 2604 737 
155777 | 128 3 | 128 2954 540 
2492417 | 2048 3 256 805 647 910 
1249 | 32 7 359713 32 | 41 32 323 896 911 
799361 | 128 3 64 514658 695 
799 361 128 3 128 505781 672 
1289! 8 | 6 371233 | 32] 5 | 32] 307451 | 792 | 
1297 | 16 | 10 83009 | 64| 31 32 79164 | 216 | 
83009 64 3 64 82 306 355 | 
1321 8 13 253633 64 11 32 72820 371 | 
1361 16 3 130657 32 5 32 58672 114 
696833 | 512 5 64 | 276791 574 
1409 | 128 3 315617 32 3 32 91 464 72 
450881 64 6 64 197 893 1254 
1623169 | 128 13 | 128] 1335737 327 
8296193 | 256 ihe 256 © me Paar 
8656897 | 2048 512 
1433| 8 | 3 733697 | 512| 31 32 8973 | 1001 | 
1481| 8 | 3 758273 | $12; 3 | 32| 309921 | 826 | 
1489 | 16 | 14 428833 32 5 32 | 340001 | 1474 
571777 | 128 5 64 318596 | 1300 
1553 | 16 3 49697 32 3 32 15122 881 
1987841 | 256 3 64 | 1532826 | 251 
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A previously prepared list of #=1(mod mq) was kept handy in the event 
that for some # it would develop that g|F. There were two such events in 
Table 1 (as noted). In each case the use of a later # showed both events to be 
“false-alarms”’ (in Table 2). Specifically they are 


q=17, m=32; p=5,44, (q|F), and p=8,164, (qtF); 
q=449, m=64; p=2,643,713, (g|F), and p= 4,367,873, (qtF). 


Thus in each case the hypothesis ‘‘g|g,” is disproved. 


Table 2. Supplementary run 














1 2 3 4 5 6 7 8 9 10 
a |u@| i |? Ma) | ¢ | < ) oe | tee 
17 | 16 3 6529 | 128 7 32 2742 9 q|F 
8161 32 7 32 2993 9 
10337 32 3 32 1692 9 q\F 
11969 64 3 32 10400 9 q|F 
15233 128 3 32 386 9 
20129 32 3 32 5055 9 q|F 
24481 32 11 32 1269 9 
26113 512 7 32 5700 9 
28289 128 6 32 7219 9 
29921 32 3 32 21479 9 q|F 
34273 32 5 32 3016 9 
37 537 32 13 32 26646 9 
40801 32 13 32 26242 9 
42433 64 10 32 11947 9 
17 16 3 6529 128 7 64 1471 3 
11969 64 3 64 5049 3 
15233 | 128 3 64 9208 3 q|F 
26113 | 512 7 64 9604 3 q|F 
28989 | 128 6 64 14715 3 
449 | 64 3 | 4367873 512 3 256 | 1876122 | 391 





























The entire output was too long for inclusion. For instance, for the first case 
(in Table 1) the output was the following array (explained in (5.13)): 


5441 3 17 3 
3234 9 
632 4998 3474 489 
630 4996 3472 487 
1453 2554 2907 —-2633 
14 16 2 0 


FFFFFFFF (Signal indicating q| F). 
Numer. Math. Bd. 2 26 
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The output was abridged for publication to show just the primitie roots g 
(of p), and 7 (of g) as well as z and x, which are essentially primitive m™ roots 
of unity modulo # and n" roots of unity modulo g. These might be of value 
for other calculations with f,. 


9. Probability of ‘‘false alarms’”’ 


We now consider the probability that for values of ZL, drawn at random, 
modulo q, it will happen that g|F. We easily see the answer to be 2 where 


(9.1) 2Q=1— (4—1/q)™"*. 
If g—1=mv/4, (v1), then for “‘small v” 
(9.2) Q ~1— exp — 1/2v <1 — exp — 1/2 ~ 39%. 


Thus, in the first case run, g=17, m=}32 the probability that q|F is roughly 
39% and for g=17, m=64, the probability is roughly 22%, (v=2). 

In Table 2, supplementary cases were taken for g=17, m=}32 with an average 
of 6/15 =40% where g|F, und for g=17, m= 64 there was an average of 2/6= 33% 
where qg|F. These values do not compare too badly with randomness. In the 
event g|g,, of course, randomness should disappear as we find q|F for each 
p =1(mod mq). 


10. The sixty-fourth root 


The most immediate challenge is f,. 

We can obtain a fairly easy estimate of 4, by asking how many different 
ideal classes could be of norm <N,=881,830. From the fact that 4;=1, we 
know that the non-principal prime ideals (if present) divide only rational primes 
of type p=64k +1. 

Now REvSCHLE’s table [4; p. 758] helps us in cases where 6<1000. When 

=—1 mod 64, p=127, 191, 383, (<1000) for which REUSCHLE displays al- 
gebraic integers « in f, for which N(a)=~. When P=1 mod 64, there are 6 
values < 1000: 


(10.1) P = 193, 257, 449, 577, 641, 769. 


REUSCHLE lists algebraic integers « of norm in &, equal (essentially) to 193 x 257, 
193 X641, 193 x 769, 769x449, 769x577. Then if « is the complex conjugate 
of «, a« is in f, and has the indicated norms in f,. Thus if 6<1000, andp=-+41 
mod 64 it follows that either the ideal factors of p are principal, or they belong 
to one of the 32 classes determined by divisors of 193 and their inverses as re- 
presented by (say) divisors of 257. 

Now GEORGE produced a list of all primes = +1 mod 64 showing exactly 
2180 primes =1 and 2207 primes =—1 mod 64, less than N,. Considering 
the fact that the nine primes <1000 have divisors that are principal or belong 
to only 32 possible classes, we find at most 


(40.2) 16 (2207 + 2180 — 9) + 32 = 70,048 + 32 


possible non-principal inequivalent integral ideals of prime norm < N,. 
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If we consider composite norms, we see a product of three (or more) ideal 
factors makes no contribution to the list of classes since 193?> N,. If we consider 
norms consisting of products of two primes, one is < N}=937. Now there are 
25 primes Q = +1 (mod 64) from 1000 to N,/193 and 19 primes Q’= +1 (mod 64) 
from 1000 to N,/257. It is easily seen that a complete set of all non-principal 
non-equivalent ideal classes of composite norm is given by some subset (possibly 
void) of ideals of norm 193 P, 257P, 193Q, 257Q’, numbering at most 


(10.3) 16X15 X6+16X15 x5 + 16225 + 162 19 = 13,904. 
Thus 
(10.4) hg = 1+ 70,048 + 32 + 13,904 = 83,965 = H. 


The closest smaller prime =1 mod 16 is 83,921. 

If we could improve slightly on REUSCHLE by showing the prime ideal factors 
of 193 to be principal in f,, then we could write 4g< 70,049. 

In either case we should have to deal with an enormous number of primes 
g=1 mod 16, as possible divisors of 4g. There are 1011 values of these g<H 
(and 859 values of g< 70,049). We have taken care of only 26 values of g through 
1553 so far and the experiment is scarely worth continuing on present equipment, 
with high values of #=1 (mod 16q) looming ahead. 

We merely verify that h, is prime or unity. We first note that H!<1553. 
Thus 4, can be composite only if the equation (4.1) holds for g =1 mod 16 (and 
q®|h,). We now treat equation (4.1) modulo 193 by using the imput data p= 193, 
q=3, m=64. (The only difference, now, is that 7, required in equation (3.16), 
does not exist, so we must ignore values of f/ in the output (5.13)). The output 
reveals g=5, z=125=2%, 2,=2§ =87 (mod 193), and 


(10.5) T, = (2, — 1)/(z + 1) = 86/88 = 5-® (mod 193). 


Hence Tj (or t) are not perfect cubes modulo 193 (or modulo all divisors of 193). 
The fifth (and higher) roots of t are ignored since 5*>H. 


11. Concluding Remarks 


We: must face the realization that the present type computer is just too 
slow by a factor of at least 100 for final results on f,, the first field that lies 
beyond hand calculations, and hopelessly slow for higher ones, with the present 
state of theory. | 

We still have obtained no evidence to doubt that every h,=1. 


We shall turn our attention, hopefully, to the possibility that t, be a perfect 
q‘* power in f, for some ‘‘small” g #1 (mod), in Part II of this work. 


This work was supported in part by the National Science Foundation Grant 
G-7412, and the computer time was donated by the Applied Mathematics Division 
of the Argonne National Laboratory of the Atomic Energy Commission. Special 
thanks are due to division director Dr. W. F. MILLER for his encouragement and 
support and to Mr. RoBERT LITTLE for his economical job of programming. The 
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computation was run on the GEORGE computer at Lemont, Illinois, August- 
October 1959. 


Added in proof (September 6, 1960): Recent work of the author ethablishes that t 
is never a perfect power, thus strengthening these results, (see § 4). 
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Taschenbuch fiir Rechenautomaten 


Ein ,,Taschenbuch fiir Rechenautomaten“, das in. fiinf oder mehr Banden erscheinen soll, 
befindet sich innerhalb der Reihe ,,Grundlehren der Mathematischen Wissenschaften“ (Gelbe 
Sammlung) in Vorbereitung. Herausgeber sind 


F. L. BAvER, Mainz 

A. S. HousEHOLDER, Oak Ridge 
F. W. J. OLvER, Teddington 

H. RuTIsHAvuseEr, Ziirich 

K. SAMELSON, Mainz 

R. Saver, Miinchen 

E. STIEFEL, Ziirich 


Zweck des Taschenbuches ist die Zusammenstellung erprobter Algorithmen fiir Rechnungen 
verschiedenster Art: Lésung endlicher Gleichungen und Funktionalgleichungen, numerische 
Approximationen, Berechnung spezieller Funktionen usw. Die Algorithmen werden in der 
Symbolsprache ALGOL abgefaBt und sind daher mit geeignetem Formeliibersetzer fiir alle 
Rechenmaschinen verwendbar. Auch ohne Formeliibersetzer kénnen sie als Programmie- 
rungsmodelle benutzt werden. Der Text der Bande ist englisch. 


Zur Zeit sind fiir die einzelnen Bande folgende Themen vorgesehen: Band 1a beschreibt 
die Symbolsprache ALGOL und ihren Gebrauch, Band 1b die Struktur der Formeliibersetzer. 
Diese einfiihrenden Bande enthalten damit als einzige nicht vorwiegend Algorithmen. Band 2 
wird der Lésung endlicher linearer und nichtlinearer Gleichungen einschlieBlich der Bestim- 
mung der Eigenwerte und Eigenvektoren von Matrizen gewidmet sein. Band 3 soll Funk- 
tionalgleichungen, insbesondere gewéhnliche und partielle Differentialgleichungen und Inte- 
gralgleichungen behandeln, Band 4 Numerische Approximationen und Band 5 die Berechnung 
spezieller Funktionen. Méglicherweise bleiben gewisse Algorithmen, z.B. zur Lésung von 
Ungleichungen, fiir das Programmieren, fiir statistische Berechnungen und 4hnliches einem 
sechsten Band vorbehalten. Zu jedem Algorithmus gehért ein ausfiihrlicher Begleittexé,; der 
iiber Zeitaufwand, Genauigkeit und Reichweite informiert und es erméglicht, die Brauchbar- 
keit eines Algorithmus fiir ein vorgegebenes Problem zu beurteilen. Es werden nur erprobte 
Algorithmen veréffentlicht. 

Vor Erscheinen der einzelnen Bande werden die Algorithmen jeweils in Teilbeitragen inner- 
halb der Zeitschrift ,,Numerische Mathematik‘ veréffentlicht, um sie médglichst bald der 
Allgemeinheit zuganglich zu machen und um vor der endgiiltigen Fassung des Bandes eventuell 
noch zusatzliche Informationen und Berichtigungen von Lesern zu erhalten. 


Beitrage zum Taschenbuch werden freundlichst erbeten. Allerdings kénnen zunachst nur 
erprobte Algorithmen angenommen werden, denen genaue Angaben iiber Art und Umfang 
der Erprobung, Genauigkeitsschatzung und ein Erfahrungsbericht beigefiigt sind. Nicht 
erprobte Algorithmen werden nicht ipso facto zuriickgewiesen, kénnen aber erst nach erfolgter 
Erprobung verwendet werden. Hinweise auf bereits veréffentlichte Algorithmen werden 
ebenfalis erbeten. Eventuelle Beitrage sind an die oben genannten Herausgeber zu richten. 
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